source: trunk/GDE/RAxML/leaveDropping.c

Last change on this file was 19480, checked in by westram, 17 months ago
File size: 40.0 KB
Line 
1#include <stdlib.h>
2#include <stdio.h>
3#include <string.h>
4#include <math.h>
5#include <assert.h>
6#include "axml.h"
7
8
9/* the define below makes the algorithm run much faster and should therefore always be used :-) */
10
11#define _SORTED
12
13/*
14  The define below is usually not required, Andre just
15  used it, to have a one to one comparison/verification with Nick's original
16  python script.
17  The only thing this changes is that taxa are sorted in exactly the same way as
18  in Nick's script. In some cases, since we are dealing with a greedy algorithm
19  the taxon order has an impact on the dropsets size and composition.
20
21
22  #define _COMPARE_TO_NICK
23
24*/
25
26extern char run_id[128];
27extern char workdir[1024];
28extern double masterTime;
29extern const unsigned int mask32[];
30
31
32
33typedef struct hash_el
34{
35  unsigned int fullKey;
36  void *value;
37  struct hash_el *next;
38} HashElem;
39
40typedef struct hash_table
41{
42  unsigned int tableSize;
43  unsigned int entryCount;
44  void *commonAttributes;
45  unsigned int (*hashFunction)(struct hash_table *h, void *value);
46  boolean (*equalFunction)(struct hash_table *hashtable, void *entrA, void *entryB);
47  HashElem **table;
48} HashTable;
49
50typedef struct 
51{
52  HashTable *hashTable;
53  HashElem *hashElem;
54  unsigned int index;
55} HashTableIterator;
56
57
58
59typedef struct profile_elem
60{
61  /* TODO should better be named taxaSet and treeSet */
62  unsigned int *bitVector;
63  unsigned int *treeVector;
64  unsigned int treeVectorSupport;
65  unsigned int id;
66} ProfileElem;
67
68typedef struct
69{
70  unsigned int *dropset;
71  unsigned int sizeOfDropset;
72  unsigned int numberOfMergingPairs; 
73} Dropset;
74
75typedef struct 
76{
77  unsigned int dropsetVectorLength;
78  unsigned int bipartitionVetorLength;
79  unsigned int *randForTaxa;
80} DropsetAttr;
81
82typedef struct 
83{
84  unsigned int bitVectorLength; 
85  unsigned int treeVectorLength; 
86  unsigned int *randForTaxa;    /* random numbers to hash the vectors */
87  unsigned int lastByte;                /* the padding bits */
88} ProfileElemAttr;
89
90
91typedef struct 
92{
93  void *arrayTable;
94  void *commonAttributes; 
95  unsigned int length;
96} Array;
97
98
99#define FLIP_NTH_BIT(bitVector,n) (bitVector[n / MASK_LENGTH] |= mask32[ n % MASK_LENGTH ])
100#define UNFLIP_NTH_BIT(bitVeoctor,n) (bitVector[n / MASK_LENGTH] &= ~mask32[ n % MASK_LENGTH ])
101#define NTH_BIT_IS_SET(bitVector,n) (bitVector[n / MASK_LENGTH] & mask32[n % MASK_LENGTH])
102
103
104typedef struct 
105{
106  List *elems;
107  unsigned int length;
108} SparseBitVector;
109
110/* SWITCHES */
111
112#define PROFILE_TABLE(x,y) createHashTable(x,y,randomProfileElemHashFunction, bitVectorEqual)
113#define DROPSET_TABLE(x,y) createHashTable(x,y,randomDropsetHashFunctionSparse, dropsetEqualFunctionSparse)
114#define SORT_ARRAY(array, elemType, sortFunction) qsort(array->arrayTable, array->length, sizeof(elemType*), sortFunction)
115#define HASH_TABLE_SIZE(x) (x * LOG(x))
116
117/* END */
118
119unsigned int genericBitCount(unsigned int* bitVector, unsigned int bitVectorLength)
120{
121  unsigned int 
122    i, 
123    result = 0;
124
125  for(i = 0; i < bitVectorLength; i++)
126    result += BIT_COUNT(bitVector[i]);
127 
128  return result; 
129}
130
131
132#ifdef _ANDRE_UNUSED_FUNCTIONS
133
134/* stuff added by Andre that is not used however */
135
136static unsigned int naiveHashFunctionAdapterVersion(unsigned int *bv, unsigned int length)
137{
138  unsigned int 
139    i, 
140    result = 0; 
141 
142  for( i = 0; i < length; ++i)
143    result ^= bv[i];
144 
145  /* return enhanceHash(result); */
146 
147  return result; 
148}
149
150
151
152static void *searchHashTableSeq(HashTable *hashtable, void *value, unsigned int hashValue)
153{
154  unsigned int 
155    position = hashValue % hashtable->tableSize;
156 
157  HashElem
158    *elem, 
159    *result = (HashElem*)NULL; 
160
161  for(elem = hashtable->table[position]; 
162      elem; 
163      elem = elem->next)
164    {
165      if(elem->fullKey  == hashValue && 
166         hashtable->equalFunction(hashtable, elem->value, value))
167        {
168          result = elem->value; 
169          break;
170        }
171    }
172
173  return result;
174}
175
176static Array* dropsetHashToArray(HashTable *dropsets)
177{
178  unsigned int 
179    count = 0; 
180 
181  HashTableIterator
182    *hashTableIterator = createHashTableIterator(dropsets);
183 
184  Array
185    *array = (Array*)rax_calloc(1, sizeof(Array));
186 
187  array->arrayTable = (void *)rax_calloc(dropsets->entryCount, sizeof(Dropset*));
188  array->length = dropsets->entryCount;
189 
190  do
191    {     
192      ((Dropset**)array->arrayTable)[count++] = getCurrentValueFromHashTableIterator(hashTableIterator);
193    } 
194  while(hashTableIteratorNext(hashTableIterator));
195 
196  assert(count == dropsets->entryCount);
197 
198  rax_free(hashTableIterator);
199 
200  return array;
201}
202
203static void printDropset(unsigned int *bitVector, unsigned int numberOfTaxa, char **nameList)
204{
205  unsigned int 
206    i;
207 
208  printBothOpen("DROPSET: ");
209 
210  for(i = 0; i < numberOfTaxa; ++i)
211    if(NTH_BIT_IS_SET(bitVector, i))
212      printBothOpen("%s,", nameList[i+1]);
213 
214  printBothOpen("\t");
215}
216
217static unsigned int naiveDropsetHashFunction(HashTable *hashTable, void *value)
218{
219  unsigned int 
220    i, 
221    result = 0,
222    length = ((DropsetAttr*)hashTable->commonAttributes)->dropsetVectorLength;
223 
224  Dropset
225    *dropset = (Dropset*)value;
226 
227  for(i = 0; i < length; ++i )
228    result ^= dropset->dropset[i];
229 
230  return result;
231}
232
233static List* appendToList(void *value, List *list)
234{ 
235  List
236    *listElem = (List *)rax_calloc(1, sizeof(List));
237 
238  listElem->value = value;
239  listElem->next = list;
240 
241  return listElem;
242}
243
244static void insertIntoHashTableSeq(HashTable *hashTable, void *value, unsigned int index)
245{ 
246  HashElem
247    *hashElem = (HashElem*)rax_calloc(1, sizeof(HashElem));
248 
249  hashElem->fullKey = index;
250 
251  index =  hashElem->fullKey % hashTable->tableSize;
252 
253  hashElem->value = value;
254  hashElem->next = hashTable->table[index];
255  hashTable->table[index] = hashElem;
256}
257#endif
258
259static HashTable *createHashTable(unsigned int size, 
260                                  void *commonAttr, 
261                                  unsigned int (*hashFunction)(HashTable *hash_table, void *value),
262                                  boolean (*equalFunction)(HashTable *hash_table, void *entryA, void *entryB))
263{ 
264  static const unsigned int 
265    initTable[] = {64, 128, 256, 512, 1024, 2048, 4096, 
266                   8192, 16384, 32768, 65536, 131072, 
267                   262144, 524288, 1048576, 2097152,
268                   4194304, 8388608, 16777216, 33554432, 
269                   67108864, 134217728, 268435456, 
270                   536870912, 1073741824, 2147483648U};
271 
272  HashTable
273    *hashTable = (HashTable*)rax_calloc(1, sizeof(HashTable));
274 
275  unsigned int
276    tableSize,
277    i,
278    primeTableLength = sizeof(initTable)/sizeof(initTable[0]),
279    maxSize = (unsigned int)-1;   
280 
281  hashTable->hashFunction = hashFunction;
282  hashTable->equalFunction = equalFunction;
283  hashTable->commonAttributes = commonAttr;
284 
285  assert(size <= maxSize); 
286  for(i = 0; initTable[i] < size && i < primeTableLength; ++i);
287  assert(i < primeTableLength);
288
289  tableSize = initTable[i];
290
291  hashTable->table = (HashElem**)rax_calloc(tableSize, sizeof(HashElem*));
292  hashTable->tableSize = tableSize; 
293  hashTable->entryCount = 0; 
294
295  return hashTable;
296}
297
298
299/* NOTE: computing hashvalue outside...cannot afford to compute it
300   twice for bit vectors */
301
302static void *searchHashTable(HashTable *Hashtable, void *value, unsigned int hashValue)
303{
304  unsigned int 
305    position = hashValue % Hashtable->tableSize;
306 
307  HashElem
308    *elem, 
309    *result = (HashElem*)NULL; 
310
311
312 
313  for(elem = Hashtable->table[position];
314      elem; 
315      elem = elem->next)
316    {
317      if(elem->fullKey  == hashValue && 
318         Hashtable->equalFunction(Hashtable, elem->value, value))
319        {
320          result = elem->value; 
321          break;
322        }
323    }
324 
325
326
327  return result;
328}
329
330
331
332
333
334
335
336static void insertIntoHashTable(HashTable *hashTable, void *value, unsigned int index)
337{ 
338  HashElem
339    *hashElem = (HashElem*)rax_calloc(1, sizeof(HashElem));
340 
341  hashElem->fullKey = index;
342 
343  index =  hashElem->fullKey % hashTable->tableSize;
344   
345  hashElem->value = value;
346  hashElem->next = hashTable->table[index];
347  hashTable->table[index] = hashElem;
348}
349
350
351static void destroyHashTable(HashTable *hashTable, void (*rax_freeValue)(void *value))
352{
353  unsigned 
354    int i; 
355 
356  HashElem
357    *elemA, 
358    *elemB,
359    **table = hashTable->table;
360 
361  for(i = 0; i < hashTable->tableSize; ++i)
362    {
363      elemA = table[i];
364      while(elemA != NULL)
365        {
366          elemB = elemA; 
367          elemA = elemA->next; 
368          if(rax_freeValue)
369            rax_freeValue(elemB->value);
370          rax_free(elemB);
371        }
372
373
374     
375    }
376
377  rax_free(hashTable->commonAttributes); 
378  rax_free(hashTable->table);
379  rax_free(hashTable);
380}
381
382
383static void updateEntryCount(HashTable *hashTable)
384{
385  unsigned int 
386    i, 
387    result = 0;
388
389  for(i = 0; i < hashTable->tableSize; ++i)
390    {
391      HashElem
392        *elem = ((HashElem**)hashTable->table)[i];
393     
394      while(elem)
395        {
396          result++;
397          elem = elem->next;
398        }
399    }
400
401  hashTable->entryCount = result;
402}
403
404
405static HashTableIterator *createHashTableIterator(HashTable *hashTable)
406{
407  unsigned 
408    int i; 
409 
410  HashTableIterator
411    *hashTableIterator = (HashTableIterator*)rax_calloc(1, sizeof(HashTableIterator));
412 
413  hashTableIterator->hashTable = hashTable;
414  hashTableIterator->hashElem = NULL;
415  hashTableIterator->index = hashTable->tableSize;
416 
417  if(!hashTable->entryCount)
418    return hashTableIterator;
419
420  for(i = 0; i < hashTable->tableSize; ++i)
421    {
422      if(hashTable->table[i])
423        {
424          hashTableIterator->hashElem = hashTable->table[i];
425          hashTableIterator->index = i;
426          break;
427        }
428    }
429 
430  return hashTableIterator;
431}
432
433static boolean hashTableIteratorNext(HashTableIterator *hashTableIterator)
434{
435  unsigned int 
436    i, 
437    tableSize = hashTableIterator->hashTable->tableSize;
438 
439  HashElem
440    *next = hashTableIterator->hashElem->next;
441 
442  if(next)
443    {
444      hashTableIterator->hashElem = next;
445      return TRUE;
446    }
447 
448  /* TODO test case for this! */
449  /* TODO looks like it should be optimised ... =/ */
450 
451  i = hashTableIterator->index + 1;
452 
453  if(i >= tableSize)
454    {
455      hashTableIterator->index = i;
456      return FALSE;
457    }
458 
459  while((next = hashTableIterator->hashTable->table[i]) == ((HashElem*)NULL))
460    {
461      if( ++i >= tableSize) 
462        {
463          hashTableIterator->index = i; 
464          return FALSE;
465        }
466    }
467 
468  hashTableIterator->index = i;
469  hashTableIterator->hashElem = next;
470
471  return (next != ((HashElem*)NULL));
472}
473
474
475/* TODO what about performance of this? */
476
477static void *getCurrentValueFromHashTableIterator(HashTableIterator *hashTableIterator)
478{ 
479  return ((hashTableIterator->hashElem) 
480          ?  hashTableIterator->hashElem->value
481          : NULL);
482}
483
484/* TODO implement remove */
485/* TODO hash table change */
486
487
488static unsigned int randomProfileElemHashFunction(HashTable *hashTable, void *value)
489{
490  unsigned int 
491    i, 
492    result = 0,
493    length = ((ProfileElemAttr*)hashTable->commonAttributes)->bitVectorLength;
494     
495  ProfileElem
496    *profileElem = (ProfileElem*)value;
497
498  for(i = 0; i < length * MASK_LENGTH; ++i) 
499    if(NTH_BIT_IS_SET(profileElem->bitVector, i))
500      result ^= ((ProfileElemAttr*)hashTable->commonAttributes)->randForTaxa[i];
501 
502  return result; 
503}
504
505
506static boolean bitVectorEqual(HashTable *Hashtable, void *entryA, void *entryB)
507{
508  unsigned int 
509    i, 
510    bitVectorLength = ((ProfileElemAttr*)Hashtable->commonAttributes)->bitVectorLength;
511 
512  unsigned int
513    *a = ((ProfileElem*)entryA)->bitVector,
514    *b = ((ProfileElem*)entryB)->bitVector;
515 
516  for(i = 0; i < bitVectorLength; ++i)
517    {
518      if(a[i] != b[i])
519        return FALSE;
520    }
521 
522  return TRUE;
523}
524
525
526
527
528static Array* profileToArray(HashTable *profile, boolean updateFrequencyCount, boolean assignIds)
529{
530  HashTableIterator* 
531    hashTableIterator = createHashTableIterator(profile);
532 
533  Array
534    *result = (Array*)rax_calloc(1, sizeof(Array));
535 
536  unsigned int 
537    count = 0;
538
539  /* remember to always copy s.t. rax_free() runs w/o problems */
540 
541  ProfileElemAttr
542    *profileElemAttr;
543 
544  result->commonAttributes = (void*)rax_calloc(1, sizeof(ProfileElemAttr));
545  result->commonAttributes = memcpy(result->commonAttributes, profile->commonAttributes, sizeof(ProfileElemAttr));
546  profileElemAttr = result->commonAttributes;
547
548  result->length = profile->entryCount;
549  result->arrayTable = (void*)rax_calloc(profile->entryCount, sizeof(ProfileElem*));
550
551  if(!hashTableIterator)
552    return result;
553 
554  do
555    {
556      ProfileElem *profileElem = getCurrentValueFromHashTableIterator(hashTableIterator);
557
558      if(updateFrequencyCount)
559        profileElem->treeVectorSupport = genericBitCount(profileElem->treeVector, profileElemAttr->treeVectorLength);
560
561      if(assignIds)
562        profileElem->id = count;
563     
564      ((ProfileElem**)result->arrayTable)[count] = profileElem;
565      assert(profileElem->bitVector && profileElem->treeVector);
566      count++;
567    }
568  while(hashTableIteratorNext(hashTableIterator));
569 
570  assert(count == profile->entryCount);
571 
572  rax_free(hashTableIterator);
573 
574  return result;
575}
576
577
578/* THIS IS AN ADAPTER */
579
580static Array *convertHashtableToArray(hashtable *oldHashtable, unsigned int bitVectorLength, unsigned int treeVectorLength, 
581                                      unsigned int *randForTaxa, unsigned int lastByte)
582{ 
583  unsigned int
584    count = 0,
585    i;
586 
587  entry
588    *ent, 
589    *entTmp; 
590 
591  Array
592    *result; 
593 
594  ProfileElemAttr
595    *attr = (ProfileElemAttr*)rax_calloc(1, sizeof(ProfileElemAttr));
596 
597  attr->bitVectorLength = bitVectorLength; 
598  attr->treeVectorLength = treeVectorLength; 
599  attr->randForTaxa = randForTaxa;
600  attr->lastByte = lastByte;
601 
602  HashTable
603    *hashTable = PROFILE_TABLE(HASH_TABLE_SIZE(oldHashtable->entryCount), attr);
604 
605  for(i = 0; i < oldHashtable->tableSize; ++i)
606    {
607      ent = oldHashtable->table[i];
608      while(ent)
609        {
610          ProfileElem
611            *profileElem = (ProfileElem*)rax_calloc(1, sizeof(ProfileElem));
612         
613          profileElem->bitVector = ent->bitVector;
614          profileElem->treeVector = ent->treeVector;
615
616          insertIntoHashTable(hashTable, profileElem, 
617                              hashTable->hashFunction(hashTable, profileElem));
618         
619          entTmp = ent->next;
620          rax_free(ent);
621          ent = entTmp;
622        }
623    }
624 
625  rax_free(oldHashtable);
626
627  updateEntryCount(hashTable);
628   
629 
630  for(i = 0; i < hashTable->tableSize; ++i)
631    {   
632      HashElem
633        *elem = ((HashElem**)hashTable->table)[i];
634     
635      while(elem)
636        {
637          ((ProfileElem*)elem->value)->id = count++;
638          elem = elem->next;
639        }
640    }
641 
642  assert(count == hashTable->entryCount);
643
644  result = profileToArray(hashTable, TRUE, TRUE);
645  destroyHashTable(hashTable, NULL);
646 
647  return result;
648}
649
650
651#ifdef _SORTED
652
653static int sortProfileElems(const void *a, const void *b)
654{
655  unsigned int aFreq = (*(ProfileElem**)a)->treeVectorSupport,
656    bFreq = (*(ProfileElem**)b)->treeVectorSupport;
657 
658  if(aFreq < bFreq)
659    return 1;
660  else if(bFreq < aFreq)
661    return -1;
662  else
663    return 0;
664}
665
666#endif
667
668
669static Array *getInfrequentBipartitions(Array *oldArray, unsigned int threshold) 
670{
671  Array
672    *array = (Array *)rax_calloc(1, sizeof(Array)); 
673 
674  ProfileElem** 
675    oldArrayTable = oldArray->arrayTable;
676
677  unsigned int 
678    i, 
679    numberInfrequent = 0, 
680    count = 0; 
681
682  for(i = 0; i < oldArray->length; ++i)
683    numberInfrequent += (oldArrayTable[i]->treeVectorSupport <= threshold);
684
685  array->commonAttributes = (void*)rax_calloc(1, sizeof(ProfileElemAttr));
686  memcpy(array->commonAttributes, oldArray->commonAttributes, sizeof(ProfileElemAttr));
687  array->length = numberInfrequent;
688  array->arrayTable = (void*)rax_calloc(numberInfrequent, sizeof(ProfileElem*));
689 
690  for(i = 0; i < oldArray->length; ++i)
691    {
692      ProfileElem *elem = ((ProfileElem**)oldArray->arrayTable)[i];
693      assert(elem->bitVector && elem->treeVector); 
694      if( elem->treeVectorSupport <= threshold)
695        ((ProfileElem**)array->arrayTable)[count++] = elem; 
696    }
697
698  assert(count == numberInfrequent); 
699 
700#ifdef _SORTED
701  SORT_ARRAY(array, ProfileElem, sortProfileElems);
702#endif
703
704  return array;
705}
706
707
708static boolean isUnionOfTreesAboveThreshold( const ProfileElemAttr *profileElemAttr, 
709                                             const ProfileElem *elemA, const ProfileElem *elemB, 
710                                             unsigned int frequencyThreshold)
711{ 
712  unsigned int 
713    i, 
714    count = 0,
715    length = profileElemAttr->treeVectorLength; 
716
717  for(i = 0; i < length; ++i )
718    count += BIT_COUNT(elemA->treeVector[i] | elemB->treeVector[i]);
719 
720  return (count > frequencyThreshold);
721}
722
723
724
725
726
727static Dropset *getBestDropset(HashTable *dropsets
728
729#ifdef _COMPARE_TO_NICK
730                               ,char **nameList, unsigned int numberOfTaxa
731#endif
732                               
733                               )
734{
735 
736#ifdef _COMPARE_TO_NICK
737  unsigned int j;   
738#endif
739 
740  int 
741    bestImpact = 0;
742 
743  unsigned int 
744    i; 
745 
746  Dropset
747    *bestDropset = (Dropset *)NULL;
748 
749  HashTableIterator
750    *hashTableIterator = createHashTableIterator(dropsets);
751 
752  if(!getCurrentValueFromHashTableIterator(hashTableIterator))
753    return NULL;
754
755  do
756    {
757      Dropset
758        *dropset = getCurrentValueFromHashTableIterator(hashTableIterator);
759     
760      unsigned int 
761        droppedTaxa = dropset->sizeOfDropset,
762        gainedBips = dropset->numberOfMergingPairs;
763     
764      int 
765        impact = gainedBips - droppedTaxa;
766     
767      if(impact > bestImpact)
768        {
769          printBothOpen("found better dropset [gained-dropped=impact] %d\t-\t%d\t=%d\n", gainedBips, droppedTaxa, impact);
770          bestImpact = impact;
771          bestDropset = dropset;
772        }
773     
774#ifdef _COMPARE_TO_NICK
775      else if(impact > 0 && impact == bestImpact)
776        {         
777          /* resolve by number of taxa */
778          unsigned int 
779            fromBest = bestDropset->sizeOfDropset,
780            fromCurrent = dropset->sizeOfDropset;
781         
782          if(fromCurrent < fromBest)
783            {
784              printBothOpen("[SIZE] ");
785              printDropset(dropset->dropset, numberOfTaxa, nameList);
786              printBothOpen(" is shorter than ");
787              printDropset(bestDropset->dropset, numberOfTaxa, nameList);
788              bestDropset = dropset;
789              printBothOpen("\n");
790            }
791          else if(fromBest == fromCurrent)
792            {
793              for(j = 0; j < fromCurrent; ++j)
794                {                 
795                  if( bestDropset->dropset[j] < dropset->dropset[j] )
796                    {
797                      printBothOpen("[LEX] ");
798                      printDropset(bestDropset->dropset, numberOfTaxa, nameList);
799                      printBothOpen("  less than ");
800                      printDropset(dropset->dropset, numberOfTaxa, nameList);
801                      printBothOpen("\n");
802                      break;
803                    }
804                  else if( dropset->dropset[j] < bestDropset->dropset[j] )
805                    {
806                      printBothOpen("[LEX] ");
807                      printDropset(dropset->dropset, numberOfTaxa, nameList);
808                      printBothOpen("  less than ");
809                      printDropset(bestDropset->dropset, numberOfTaxa, nameList);
810                      bestDropset = dropset;
811                      printBothOpen("\n");
812                      break;
813                    }
814                }
815            }
816        }
817#endif
818    } 
819  while(hashTableIteratorNext(hashTableIterator));
820
821  rax_free(hashTableIterator); 
822 
823  /* we did not find anything */
824  if(!bestDropset)
825    return bestDropset;
826
827  /* else convert to non-sparse representation */
828  unsigned int 
829    *bvPtr = bestDropset->dropset; 
830 
831  bestDropset->dropset = (unsigned int*)rax_calloc(((DropsetAttr*)dropsets->commonAttributes)->dropsetVectorLength, sizeof(unsigned int));
832
833  for(i = 0; i < bestDropset->sizeOfDropset; ++i)
834    FLIP_NTH_BIT(bestDropset->dropset, bvPtr[i]);
835 
836  rax_free(bvPtr);
837 
838  return bestDropset;
839}
840 
841 
842static void insertBipartitionPairDropset_helper(HashTable *dropsets,
843                                                unsigned int *diff, unsigned int diffCount)
844{
845  Dropset
846    *dropset = (Dropset*)rax_calloc(1, sizeof(Dropset)),
847    *foundInHashTable;
848  unsigned int 
849    i, 
850    ctr = 0;
851 
852  dropset->dropset = (unsigned int*)rax_calloc(diffCount, sizeof(unsigned int));
853 
854  for(i = 0; i < ((DropsetAttr*)dropsets->commonAttributes)->dropsetVectorLength * MASK_LENGTH; ++i)
855    {
856      if(NTH_BIT_IS_SET(diff, i))
857        {
858          dropset->dropset[ctr] = i;
859          ctr++;
860          if(ctr == diffCount)
861            break;
862        }
863    }
864 
865  rax_free(diff);
866
867  assert(ctr == diffCount);
868
869  dropset->numberOfMergingPairs = 1;
870  dropset->sizeOfDropset = diffCount;
871  unsigned int hashValue = dropsets->hashFunction(dropsets, dropset); 
872 
873  foundInHashTable = (Dropset*)searchHashTable(dropsets, dropset, hashValue);
874 
875  if(!foundInHashTable)
876    insertIntoHashTable(dropsets, dropset, hashValue);
877  else
878    {
879      foundInHashTable->numberOfMergingPairs++;
880           
881      rax_free(dropset->dropset);
882      rax_free(dropset);
883    }
884}
885
886
887static void insertBipartitionPairDropset(HashTable *dropsets,
888                                         const ProfileElemAttr  *profileElemAttr, 
889                                         const ProfileElem *elemA, const ProfileElem *elemB,
890                                         const unsigned int *droppedTaxa )
891{
892  unsigned int 
893    i, 
894    diffCount1 = 0, 
895    diffCount2 = 0,
896    length = profileElemAttr->bitVectorLength,
897    *diff1 = (unsigned int*)rax_calloc(length, sizeof(unsigned int)), 
898    *diff2 = (unsigned int*)rax_calloc(length, sizeof(unsigned int)),
899    lastByte = profileElemAttr->lastByte;
900
901 
902  for(i = 0; i < length; ++i )
903    {
904      diff1[i] = ( elemA->bitVector[i] ^ elemB->bitVector[i] ) ;
905
906      if(i == length - 1)
907        diff2[i] = (elemA->bitVector[i] & ~ droppedTaxa[i])  ^ ~ (elemB->bitVector[i] | droppedTaxa[i] | lastByte); 
908      else
909        diff2[i] = (elemA->bitVector[i] & ~ droppedTaxa[i]) ^ ~ (elemB->bitVector[i] | droppedTaxa[i]); 
910    }
911
912 
913  diffCount1 = genericBitCount(diff1,length);
914  diffCount2 = genericBitCount(diff2,length);
915 
916  if (!(diffCount1 && diffCount2)) 
917    printf("problem with bip %d and %d\n", elemA->id, elemB->id);
918 
919  assert( diffCount1 && diffCount2 ) ; 
920
921  if(diffCount1 < diffCount2)
922    {
923      insertBipartitionPairDropset_helper(dropsets, diff1, diffCount1);
924      rax_free(diff2);
925    }
926  else if(diffCount1 > diffCount2)
927    {
928      insertBipartitionPairDropset_helper(dropsets, diff2, diffCount2);
929      rax_free(diff1);
930    }
931  else
932    {
933      insertBipartitionPairDropset_helper(dropsets, diff1, diffCount1);
934      insertBipartitionPairDropset_helper(dropsets, diff2, diffCount2);
935    }
936}
937
938
939
940
941static unsigned int randomDropsetHashFunctionSparse(HashTable *hashTable, void *value)
942{
943  Dropset
944    *dropset = (Dropset*)value;
945 
946  unsigned int 
947    i, 
948    result = 0,
949    length = dropset->sizeOfDropset,
950    *randForTaxa = ((DropsetAttr*)hashTable->commonAttributes)->randForTaxa;
951
952  for(i = 0; i < length; ++i)
953    result ^= randForTaxa[dropset->dropset[i]];
954
955  return result; 
956}
957
958
959static boolean dropsetEqualFunctionSparse(HashTable *hashTable, void *entryA, void *entryB)
960{
961    (void)hashTable; // suppress warning
962  unsigned int i,
963    aLength = ((Dropset*)entryA)->sizeOfDropset,
964    bLength = ((Dropset*)entryB)->sizeOfDropset;
965 
966  unsigned int
967    *a = ((Dropset*)entryA)->dropset,
968    *b = ((Dropset*)entryB)->dropset;
969 
970  if(aLength != bLength)
971    return FALSE;
972 
973  for(i = 0; i < aLength; i++)
974    if(a[i] != b[i])
975      return FALSE;
976
977  return TRUE;
978}
979
980
981static HashTable* potentialProfileDropsets(Array *infrequentBipartitions,               
982                                           unsigned int frequencyThreshold,
983                                           unsigned int *droppedTaxa )
984{
985  DropsetAttr
986    *dropsetAttr = (DropsetAttr*)rax_calloc(1, sizeof(DropsetAttr));
987 
988  dropsetAttr->dropsetVectorLength = ((ProfileElemAttr*)infrequentBipartitions->commonAttributes)->bitVectorLength;
989  dropsetAttr->bipartitionVetorLength = infrequentBipartitions->length;
990  dropsetAttr->randForTaxa = ((ProfileElemAttr*)infrequentBipartitions->commonAttributes)->randForTaxa;
991
992  HashTable
993    *dropsets = DROPSET_TABLE(HASH_TABLE_SIZE(infrequentBipartitions->length), dropsetAttr);
994 
995  ProfileElem
996    *elemA, 
997    *elemB;   
998 
999  unsigned int 
1000    firstCount, 
1001    secondCount; 
1002 
1003  for(firstCount = 0; firstCount < infrequentBipartitions->length; ++firstCount)
1004    {
1005      elemA = ((ProfileElem**)(infrequentBipartitions->arrayTable))[firstCount];
1006      for(secondCount = firstCount + 1; secondCount < infrequentBipartitions->length; ++secondCount)
1007        {
1008          elemB = ((ProfileElem**)(infrequentBipartitions->arrayTable))[secondCount];
1009          assert(elemB->treeVector && elemB->bitVector);
1010
1011          if( elemA->treeVectorSupport + elemB->treeVectorSupport > frequencyThreshold )
1012            {
1013              if( isUnionOfTreesAboveThreshold(infrequentBipartitions->commonAttributes, elemA, elemB, frequencyThreshold))
1014                insertBipartitionPairDropset(dropsets, infrequentBipartitions->commonAttributes, elemA, elemB, droppedTaxa);
1015            }
1016#ifdef _SORTED
1017          else
1018            break;
1019#endif
1020        }
1021    }
1022 
1023  /* due to parallelisation, we have to compute the entry count a posteriori */
1024  /* TODO rewrite this step for the sequential code */
1025
1026  updateEntryCount(dropsets);
1027
1028  return dropsets; 
1029}
1030
1031
1032static void destroyProfileElem(ProfileElem *profileElem)
1033{
1034  rax_free(profileElem->bitVector);
1035  rax_free(profileElem->treeVector);
1036  rax_free(profileElem);
1037}
1038
1039
1040static void destroyDropset(void *dropset_) 
1041{
1042  Dropset
1043    *dropset = (Dropset*)dropset_;
1044 
1045  rax_free(dropset->dropset);
1046  rax_free(dropset);
1047}
1048
1049
1050static List* getListOfConsensusBips(Array *allBipartitions, unsigned int frequencyThreshold)
1051{
1052  List* 
1053    result = (List*)NULL;       
1054 
1055  unsigned int 
1056    i;
1057 
1058  ProfileElem
1059    *elem; 
1060
1061  for(i = 0; i < allBipartitions->length; ++i)
1062    {
1063      elem = ((ProfileElem**)allBipartitions->arrayTable)[i];
1064      if(elem->treeVectorSupport > frequencyThreshold)
1065        {
1066          List
1067            *tmp = (List*)rax_calloc(1, sizeof(List));
1068         
1069          tmp->value = elem;
1070          tmp->next = result;   
1071          result = tmp;
1072        }
1073    }
1074 
1075  return result;
1076}
1077 
1078
1079static unsigned int getLengthOfList(List* list)
1080{ 
1081  unsigned int 
1082    result = 0; 
1083
1084  for( ; list ; list = list->next)
1085    result++;
1086   
1087  return result;
1088}
1089
1090
1091
1092
1093
1094static HashTable *updateAndInsertElem(ProfileElem *elem, HashTable *hashTable, unsigned int *droppedTaxa, ProfileElemAttr* profileElemAttr)
1095{
1096  unsigned int 
1097    j; 
1098 
1099  ProfileElem
1100    *foundInHashTable, *foundComplementInHashTable,
1101    *complement = (ProfileElem*)rax_calloc(1, sizeof(ProfileElem));
1102 
1103  unsigned int 
1104    lastByte = ((ProfileElemAttr*)hashTable->commonAttributes)->lastByte;
1105 
1106  complement->bitVector = (unsigned int *)rax_calloc(profileElemAttr->bitVectorLength, sizeof(unsigned int));
1107
1108  for(j = 0; j < profileElemAttr->bitVectorLength; ++j )
1109    {
1110      if(j == profileElemAttr->bitVectorLength - 1)
1111        {
1112          elem->bitVector[j] &= ~ ( droppedTaxa[j] | lastByte );
1113          complement->bitVector[j] = ~ (elem->bitVector[j] | droppedTaxa[j] | lastByte);
1114        }
1115      else
1116        {
1117          elem->bitVector[j] &= ~ droppedTaxa[j];
1118          complement->bitVector[j] = ~ (elem->bitVector[j] | droppedTaxa[j] );
1119        }
1120    }
1121 
1122  /* bipartition vanishes */
1123 
1124  unsigned int 
1125    numberOfElements = genericBitCount(elem->bitVector, profileElemAttr->bitVectorLength),
1126    numberElementsInComplement = genericBitCount(complement->bitVector, profileElemAttr->bitVectorLength); 
1127 
1128  if(numberOfElements < 2 || numberElementsInComplement < 2)
1129    {
1130      destroyProfileElem(elem);     
1131      rax_free(complement->bitVector);
1132      rax_free(complement);
1133      return hashTable;
1134    }
1135 
1136  unsigned int 
1137    hashValue = hashTable->hashFunction(hashTable, elem);
1138
1139  foundInHashTable = searchHashTable(hashTable, elem, hashValue);
1140 
1141  if(!foundInHashTable)
1142    {
1143      unsigned int 
1144        complementHashValue = hashTable->hashFunction(hashTable, complement);
1145     
1146      foundComplementInHashTable = searchHashTable(hashTable, complement, complementHashValue);
1147     
1148      if(foundComplementInHashTable)
1149        foundInHashTable = foundComplementInHashTable;
1150    }
1151
1152  if(foundInHashTable)
1153    {     
1154      for (j = 0; j < profileElemAttr->treeVectorLength; ++j )
1155        foundInHashTable->treeVector[j] |= elem->treeVector[j];     
1156     
1157      foundInHashTable->treeVectorSupport = genericBitCount(foundInHashTable->treeVector, profileElemAttr->treeVectorLength);
1158     
1159      destroyProfileElem(elem);
1160      rax_free(complement->bitVector);
1161      rax_free(complement);
1162    }
1163  else
1164    {
1165      insertIntoHashTable(hashTable, elem, hashValue);
1166      rax_free(complement->bitVector);
1167      rax_free(complement);
1168    } 
1169
1170  return hashTable;
1171}
1172
1173
1174static Array* restrictProfile( Array *infrequentBipartitions, List *consensusBipartitions, unsigned int *droppedTaxa )
1175{
1176  unsigned int 
1177    i; 
1178 
1179  ProfileElemAttr
1180    *profileElemAttr = (ProfileElemAttr *)rax_calloc(1, sizeof(ProfileElemAttr));
1181 
1182  profileElemAttr = memcpy(profileElemAttr, infrequentBipartitions->commonAttributes, sizeof(ProfileElemAttr));
1183 
1184  HashTable
1185    *hashTable = PROFILE_TABLE(HASH_TABLE_SIZE(infrequentBipartitions->length + getLengthOfList(consensusBipartitions)),
1186                               profileElemAttr);
1187 
1188  unsigned int 
1189    lengthOfConsensus = getLengthOfList(consensusBipartitions);
1190 
1191  Array
1192    *result; 
1193 
1194  List
1195    *listElem, 
1196    *listTmp;
1197
1198  Array
1199    *tmpArray = (Array *)rax_calloc(1, sizeof(Array));
1200
1201 
1202  tmpArray->arrayTable = (void *)rax_calloc(infrequentBipartitions->length + lengthOfConsensus, sizeof(ProfileElem*));
1203  tmpArray->commonAttributes = (void *)rax_calloc(1, sizeof(ProfileElemAttr));
1204  tmpArray->commonAttributes = memcpy(tmpArray->commonAttributes, infrequentBipartitions->commonAttributes, sizeof(ProfileElemAttr));
1205
1206  for(i = 0; i < infrequentBipartitions->length; ++i)
1207    ((ProfileElem**)tmpArray->arrayTable)[i] = ((ProfileElem**)infrequentBipartitions->arrayTable)[i];
1208
1209  tmpArray->length = infrequentBipartitions->length;
1210
1211  rax_free(infrequentBipartitions->commonAttributes);
1212  rax_free((ProfileElem**)infrequentBipartitions->arrayTable);
1213  rax_free(infrequentBipartitions);
1214
1215  listElem = consensusBipartitions;
1216 
1217  while(listElem)
1218    {
1219      ((ProfileElem**)tmpArray->arrayTable)[tmpArray->length] = (ProfileElem*)listElem->value;
1220      tmpArray->length++;
1221      listTmp = listElem;
1222      listElem = listElem->next;
1223      rax_free(listTmp);
1224    } 
1225 
1226  for(i = 0; i < tmpArray->length; ++i)
1227    {
1228      hashTable = updateAndInsertElem(((ProfileElem**)tmpArray->arrayTable)[i],
1229                                      hashTable,
1230                                      droppedTaxa,
1231                                      profileElemAttr);
1232    }
1233
1234  rax_free(tmpArray->arrayTable);
1235  rax_free(tmpArray->commonAttributes);
1236  rax_free(tmpArray);
1237
1238  updateEntryCount(hashTable);
1239
1240  result = profileToArray(hashTable, TRUE, FALSE);
1241
1242  destroyHashTable(hashTable, NULL);
1243  return result;
1244}
1245
1246
1247/**
1248   Computes in a greedy manner a set of rogue taxa.
1249
1250  @return the dropset
1251  @param  profile            -- the hashtable used in bipartitionList.c
1252  @param  tree               -- the standard tree
1253  @param resultBipartitions  -- bipartitions after dropping the rogue
1254                                taxa. necessary part of the output!
1255 */
1256
1257static unsigned int* determineGreedyDropset(hashtable *profile, tree *Tree, Array **resultBipartitions)
1258{   
1259
1260  unsigned int 
1261    numberOfTaxa = Tree->mxtips,
1262    numberOfTrees = Tree->numberOfTrees,
1263    frequencyThreshold,   
1264    i,
1265    treeVectorLength = GET_BITVECTOR_LENGTH(numberOfTrees),
1266    bitVectorLength = GET_BITVECTOR_LENGTH(numberOfTaxa),
1267    *droppedTaxa = (unsigned int*)rax_calloc(bitVectorLength, sizeof(unsigned int)); 
1268 
1269  boolean
1270    droppedTaxaThisRound = FALSE;
1271 
1272  HashTable
1273    *dropsets; 
1274 
1275  Dropset
1276    *bestDropset; 
1277 
1278  Array
1279    *allBipartitions;
1280
1281 
1282
1283  /* initialise random numbers: one for each taxon  */
1284  unsigned int 
1285    *randForTaxa = (unsigned int*)rax_calloc(numberOfTaxa, sizeof(unsigned int));
1286 
1287  for(i = 0; i < numberOfTaxa; ++i) 
1288    randForTaxa[i] = rand();
1289 
1290  switch(Tree->consensusType)
1291    {           
1292    case MR_CONSENSUS: 
1293      frequencyThreshold =  numberOfTrees / 2;
1294      break;
1295    case STRICT_CONSENSUS:
1296      frequencyThreshold =  numberOfTrees - 1;
1297      break;
1298    case MRE_CONSENSUS:
1299    default:   
1300      assert(0);
1301    }
1302
1303  /*
1304     prepare a compensator for the padding bits (bits that are
1305     necessarily in the bitvector but do not represent taxa) 
1306  */
1307 
1308  unsigned int 
1309    lastByte = 0;
1310 
1311  for(i = numberOfTaxa; i < MASK_LENGTH * bitVectorLength; ++i)
1312    lastByte |= mask32[i % MASK_LENGTH];
1313
1314  assert(numberOfTrees > 0 && numberOfTaxa > 0 &&  bitVectorLength > 0 && treeVectorLength > 0 && frequencyThreshold > 0);
1315 
1316  /*
1317     caling adapter: adapts the old hashtable to the new implementation
1318  */ 
1319 
1320  allBipartitions = convertHashtableToArray(profile, bitVectorLength, treeVectorLength, randForTaxa, lastByte);
1321
1322  printBothOpen("Tree->consensusType=%d, numberOfTaxa=%d, bitVectorLength=%d, numberOfTrees=%d, treeVectorLength=%d, frequencyThreshold=%d, profile->entryCount=%d\n\n",
1323                Tree->consensusType, numberOfTaxa, bitVectorLength,
1324                numberOfTrees, treeVectorLength, frequencyThreshold,allBipartitions->length); 
1325
1326 
1327
1328  do 
1329    {
1330      printBothOpen("================================================================\n");
1331
1332     
1333      Array
1334        *infrequentBipartitions = getInfrequentBipartitions(allBipartitions, frequencyThreshold);
1335         
1336      List
1337        *consensusBipartitions = getListOfConsensusBips(allBipartitions, frequencyThreshold);
1338     
1339      assert(allBipartitions->length == infrequentBipartitions->length + getLengthOfList(consensusBipartitions));       
1340   
1341      printBothOpen("divided bips: %d = %d infreq + %d consensus\n", allBipartitions->length, infrequentBipartitions->length, getLengthOfList(consensusBipartitions));
1342
1343      rax_free((ProfileElem**)allBipartitions->arrayTable);
1344      rax_free(allBipartitions->commonAttributes);
1345      rax_free(allBipartitions);
1346     
1347      /*
1348         NOTE: this below would be possible and might make sense in some
1349         situations. But as there is no equivalent in Nick's script and
1350         this somehow contradicts the philosophy stated in the paper (to
1351         remove all the rogues), I comment it out.
1352         It would of advantage, if the algorithm is used to get a better
1353         resolved tree. However, when we stop via this criteria, rogues
1354         might remain in the tree. If it is the target to yield a high
1355         number of branches with the highest possible support, then it
1356         would be better, if those rogues were dropped as well.
1357      */
1358     
1359      /*
1360         if(getLengthOfList(consensusBipartitions) >= numberOfTaxa - 3 - genericBitCount(droppedTaxa))
1361         {
1362         printBothOpen("we have enough bips: %d / %d\n", getLengthOfList(consensusBipartitions), numberOfTaxa-3);       
1363         break;
1364         }
1365         else
1366         printBothOpen("we have %d bips, we need %d.\n", getLengthOfList(consensusBipartitions), numberOfTaxa-3);
1367      */
1368
1369       
1370     
1371
1372     
1373      dropsets =  potentialProfileDropsets(infrequentBipartitions, frequencyThreshold, droppedTaxa);
1374   
1375     
1376
1377#ifdef _COMPARE_TO_NICK
1378      bestDropset = getBestDropset(dropsets , Tree->nameList, numberOfTaxa );
1379#else
1380      bestDropset = getBestDropset(dropsets);
1381#endif
1382   
1383      droppedTaxaThisRound = (bestDropset != NULL);     
1384     
1385   
1386      if(droppedTaxaThisRound)
1387        {               
1388          /*
1389             update dropped taxa
1390          */
1391         
1392          for(i = 0 ; i < ((DropsetAttr*)dropsets->commonAttributes)->dropsetVectorLength; ++i)
1393            droppedTaxa[i] |= bestDropset->dropset[i];
1394         
1395          /*
1396             update profile array
1397          */
1398         
1399          allBipartitions = restrictProfile(infrequentBipartitions, consensusBipartitions, droppedTaxa );
1400        }
1401      else
1402        {
1403          allBipartitions = restrictProfile( infrequentBipartitions, consensusBipartitions, droppedTaxa );
1404          *resultBipartitions = allBipartitions;
1405        }
1406     
1407      /* just debug */
1408     
1409      if(droppedTaxaThisRound)
1410        {
1411          printBothOpen("dropping taxa: " );
1412          for(i = 0; i < numberOfTaxa; ++i)
1413            if(NTH_BIT_IS_SET(bestDropset->dropset, i))
1414              printBothOpen("%s,", Tree->nameList[i+1]);
1415          printBothOpen("\n");
1416        }         
1417     
1418      destroyHashTable(dropsets, destroyDropset);           
1419    }     
1420  while(droppedTaxaThisRound);
1421 
1422  rax_free(randForTaxa);
1423
1424  return droppedTaxa;
1425}
1426
1427static unsigned int renormalizeBipartitions(Array *bipartitions, unsigned int *droppedTaxa, unsigned int numberOfTaxa)
1428{
1429  unsigned 
1430    i, 
1431    resultIndex = 0,
1432    bitVectorLength = ((ProfileElemAttr*)bipartitions->commonAttributes)->bitVectorLength,
1433    lastByte = 0;
1434 
1435
1436  for(i = numberOfTaxa; i < MASK_LENGTH * bitVectorLength; ++i)
1437    lastByte |= mask32[i % MASK_LENGTH];
1438 
1439  /* if the reference leave did not change, nothing needs to be done */
1440 
1441  if(!NTH_BIT_IS_SET(droppedTaxa, 0))
1442    return 0; 
1443
1444  /* find adequate new leave */
1445 
1446  for(i = 0; i < numberOfTaxa; ++i)
1447    {
1448      if(!NTH_BIT_IS_SET(droppedTaxa,i))
1449        {
1450          resultIndex = i;
1451          break;
1452        }
1453    }
1454 
1455  /* is only the case if everything was dropped  */
1456 
1457  assert(resultIndex);
1458 
1459  /* invert bit-vectors according to the reference leave */
1460 
1461  for(i = 0; i < bipartitions->length; ++i)
1462    {
1463      ProfileElem
1464        *elem =  ((ProfileElem**)bipartitions->arrayTable)[i];
1465     
1466      if(NTH_BIT_IS_SET(elem->bitVector,resultIndex))
1467        {
1468          unsigned int 
1469            j;
1470         
1471          for(j = 0; j < bitVectorLength; ++j)
1472            {         
1473              elem->bitVector[j] = 
1474                (j == bitVectorLength - 1)
1475                ? ~ (elem->bitVector[j] | droppedTaxa[j] | lastByte)
1476                : ~ (elem->bitVector[j] | droppedTaxa[j]);
1477            }
1478        }
1479    }
1480
1481  return resultIndex;
1482}
1483
1484static hashtable *reconvertHashtable(Array *bipartitionArray)
1485{
1486  unsigned int
1487    i, 
1488    count = 0; 
1489 
1490  hashtable
1491    *htable = initHashTable(bipartitionArray->length); 
1492 
1493  assert(htable->tableSize >= bipartitionArray->length );
1494 
1495  for(i = 0; i < bipartitionArray->length; ++i)
1496    {
1497      ProfileElem
1498        *elem = ((ProfileElem**)bipartitionArray->arrayTable)[i];
1499     
1500      entry
1501        *ent = initEntry();
1502     
1503      ent->bitVector = elem->bitVector;
1504      ent->treeVector = elem->treeVector;
1505      ent->bipNumber = count++;
1506      htable->table[i] = ent; 
1507      rax_free(elem);
1508    }
1509 
1510  htable->entryCount = count;
1511 
1512  return htable;
1513}
1514
1515static void pruneTaxon(tree *tr, unsigned int k)
1516{
1517
1518  assert(k > 0 && k <= ((unsigned int)(tr->mxtips)));
1519
1520  nodeptr
1521    p = tr->nodep[k],   
1522    q = p->back,
1523    q1 = q->next->back,
1524    q2 = q->next->next->back;
1525
1526
1527  hookupDefault(q1, q2, tr->numBranches);
1528
1529  tr->start = findAnyTip(q1, tr->mxtips);
1530
1531  assert(p != tr->start && q != tr->start);
1532}
1533 
1534
1535
1536
1537void computeRogueTaxa(tree *tr,  char* treeSetFileName, analdef *adef)
1538{
1539  char 
1540    dropFileName[1024];
1541
1542  hashtable
1543    *h = initHashTable(tr->mxtips * FC_INIT * 10);
1544
1545  Array
1546    **resultBipartitions = (Array **)rax_calloc(1, sizeof(Array*));
1547
1548  unsigned int 
1549    droppedTaxaNum = 0,
1550    numberOfTrees = 0, 
1551    i,   
1552    treeVectorLength, 
1553    vectorLength,
1554    **bitVectors = initBitVector(tr, &vectorLength),
1555    *droppedTaxa; 
1556
1557  FILE
1558    *treeFile = getNumberOfTrees(tr, treeSetFileName, adef),
1559    *outf;   
1560
1561  numberOfTrees = tr->numberOfTrees;     
1562
1563  assert(adef->leaveDropMode);
1564  assert(sizeof(unsigned char) == 1); 
1565 
1566  treeVectorLength = GET_BITVECTOR_LENGTH(numberOfTrees);
1567
1568  /* read the trees and process the bipartitions */
1569
1570  for(i = 1; i <= numberOfTrees; i++)
1571    {                 
1572      int 
1573        bCount = 0;
1574     
1575      treeReadLen(treeFile, tr, FALSE, FALSE, TRUE, adef, TRUE, FALSE);
1576     
1577      assert(tr->mxtips == tr->ntips);
1578     
1579      bitVectorInitravSpecial(bitVectors, tr->nodep[1]->back, tr->mxtips, vectorLength, h, (i - 1), BIPARTITIONS_BOOTSTOP, (branchInfo *)NULL,
1580                              &bCount, treeVectorLength, FALSE, FALSE);
1581     
1582      assert(bCount == tr->mxtips - 3); 
1583    } 
1584
1585
1586  rewind(treeFile); 
1587               
1588  droppedTaxa = determineGreedyDropset(h, tr, resultBipartitions);
1589       
1590
1591  renormalizeBipartitions(*resultBipartitions, droppedTaxa, tr->mxtips);
1592
1593  h = reconvertHashtable(*resultBipartitions);
1594
1595  /* TODO output for egrep */
1596   
1597  printBothOpen("> ALL dropped taxa: ");
1598 
1599  for(i = 0; i < (unsigned int)tr->mxtips; ++i)
1600    {
1601      if(NTH_BIT_IS_SET(droppedTaxa, i))
1602        {
1603          printBothOpen("%d %s,", (i+1), tr->nameList[i + 1]);
1604          droppedTaxaNum++;
1605        }
1606    }
1607  printBothOpen("\n"); 
1608   
1609
1610  printBothOpen("\nDropping %u taxa\n", droppedTaxaNum);
1611
1612  printBothOpen("\nTime for dropset calculation: %f seconds\n", gettime() - masterTime );
1613   
1614  if(droppedTaxaNum > 0)
1615    {
1616      strcpy(dropFileName, workdir);
1617      strcat(dropFileName, "RAxML_prunedTrees.");
1618      strcat(dropFileName, run_id);
1619     
1620      outf = myfopen(dropFileName, "w");
1621     
1622      for(i = 1; i <= numberOfTrees; i++)
1623        {                 
1624          unsigned int
1625            k;
1626         
1627          int
1628            tips = 0;
1629         
1630          /*printf("Tree %d\n", i);*/
1631         
1632          treeReadLen(treeFile, tr, FALSE, FALSE, TRUE, adef, TRUE, FALSE);
1633         
1634          assert(tr->mxtips == tr->ntips);
1635         
1636          for(k = 0; k < (unsigned int)tr->mxtips; k++)
1637            {
1638              if(NTH_BIT_IS_SET(droppedTaxa, k))
1639                pruneTaxon(tr, (k+1));
1640            }
1641         
1642          tips = countTips(tr->start, tr->mxtips) + countTips(tr->start->back, tr->mxtips);
1643         
1644          assert((unsigned)tips == ((unsigned)tr->mxtips - droppedTaxaNum));
1645         
1646          Tree2String(tr->tree_string, tr, tr->start->back, FALSE, TRUE, FALSE, FALSE, FALSE, adef, NO_BRANCHES, FALSE, FALSE, FALSE, FALSE);
1647          fprintf(outf, "%s", tr->tree_string);
1648         
1649          /*printf("%u %u\n", (unsigned)tips, ((unsigned)tr->mxtips - droppedTaxaNum));*/
1650        }
1651     
1652     
1653      printBothOpen("\nA tree collection, where the taxa from the dropset have been pruned in each tree\n");
1654      printBothOpen("has been written to file %s\n\n", dropFileName);
1655      fclose(outf);
1656    }
1657     
1658  printBothOpen("Total execution time: %f\n", gettime() - masterTime);
1659 
1660  fclose(treeFile); 
1661
1662  freeBitVectors(bitVectors, 2 * tr->mxtips);   
1663  rax_free(droppedTaxa);
1664  rax_free(bitVectors); 
1665  freeHashTable(h);
1666  rax_free(h);
1667
1668  exit(0);
1669}
Note: See TracBrowser for help on using the repository browser.