1 /**
2  * PLL (version 1.0.0) a software library for phylogenetic inference
3  * Copyright (C) 2013 Tomas Flouri and Alexandros Stamatakis
4  *
5  * Derived from
6  * RAxML-HPC, a program for sequential and parallel estimation of phylogenetic
7  * trees by Alexandros Stamatakis
8  *
9  * This program is free software: you can redistribute it and/or modify it
10  * under the terms of the GNU General Public License as published by the Free
11  * Software Foundation, either version 3 of the License, or (at your option)
12  * any later version.
13  *
14  * This program is distributed in the hope that it will be useful, but WITHOUT
15  * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
16  * FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License for
17  * more details.
18  *
19  * You should have received a copy of the GNU General Public License along with
20  * this program.  If not, see <http://www.gnu.org/licenses/>.
21  *
22  * For any other enquiries send an Email to Tomas Flouri
23  * Tomas.Flouri@h-its.org
24  *
25  * When publishing work that uses PLL please cite PLL
26  *
27  * @file bipartitionList.c
28  */
29 #include "mem_alloc.h"
30 
31 #ifndef WIN32
32 #include <sys/times.h>
33 #include <sys/types.h>
34 #include <sys/time.h>
35 #include <unistd.h>
36 #endif
37 
38 #include <limits.h>
39 #include <math.h>
40 #include <time.h>
41 #include <stdlib.h>
42 #include <stdio.h>
43 #include <ctype.h>
44 #include <string.h>
45 #include <stdint.h>
46 #include <assert.h>
47 
48 #include "pll.h"
49 #include "pllInternal.h"
50 
51 
52 static pllBipartitionEntry *initEntry(void);
53 static void getxnodeBips (nodeptr p);
54 static void newviewBipartitions(unsigned int **bitVectors,
55                                 nodeptr p,
56                                 int numsp,
57                                 unsigned int vectorLength,
58                                 int processID);
59 
60 static void insertHashRF(unsigned int *bitVector,
61                          pllHashTable *h,
62                          unsigned int vectorLength,
63                          int treeNumber,
64                          int treeVectorLength,
65                          hashNumberType position,
66                          int support,
67                          pllBoolean computeWRF);
68 
69 extern const unsigned int mask32[32];
70 
71 
getxnodeBips(nodeptr p)72 static void getxnodeBips (nodeptr p)
73 {
74   nodeptr  s;
75 
76   if ((s = p->next)->xBips || (s = s->next)->xBips)
77     {
78       p->xBips = s->xBips;
79       s->xBips = 0;
80     }
81 
82   assert(p->xBips);
83 }
84 
85 
initEntry(void)86 static pllBipartitionEntry *initEntry(void)
87 {
88   pllBipartitionEntry * e = (pllBipartitionEntry *)rax_malloc(sizeof(pllBipartitionEntry));
89 
90   e->bitVector     = (unsigned int*)NULL;
91   e->treeVector    = (unsigned int*)NULL;
92   e->supportVector = (int*)NULL;
93   e->bipNumber  = 0;
94   e->bipNumber2 = 0;
95   e->supportFromTreeset[0] = 0;
96   e->supportFromTreeset[1] = 0;
97   e->next       = (pllBipartitionEntry *)NULL;
98 
99   return e;
100 }
101 
cleanupHashTable(pllHashTable * h,int state)102 void cleanupHashTable(pllHashTable *h, int state)
103 {
104   unsigned int
105     k,
106     entryCount = 0,
107     removeCount = 0;
108 
109   assert(state == 1 || state == 0);
110 
111   for(k = 0, entryCount = 0; k < h->size; k++)
112     {
113       pllHashItem * start     = NULL;
114       pllHashItem * lastValid = NULL;
115 
116       pllHashItem * hitem = h->Items[k];
117       while (hitem)
118        {
119          pllBipartitionEntry *e = (pllBipartitionEntry *)(hitem->data);
120          if(state == 0)
121            {
122              e->treeVector[0] = e->treeVector[0] & 2;
123              assert(!(e->treeVector[0] & 1));
124            }
125          else
126            {
127              e->treeVector[0] = e->treeVector[0] & 1;
128              assert(!(e->treeVector[0] & 2));
129            }
130 
131          if(e->treeVector[0] != 0)
132            {
133              if(!start)
134                start = hitem;
135              lastValid = hitem;
136              hitem = hitem->next;
137            }
138          else
139            {
140              pllHashItem *tmp = hitem;
141              pllBipartitionEntry *remove = e;
142              hitem = hitem->next;
143 
144              removeCount++;
145 
146              if(lastValid) lastValid->next = hitem;
147 
148              if(remove->bitVector)     rax_free(remove->bitVector);
149              if(remove->treeVector)    rax_free(remove->treeVector);
150              if(remove->supportVector) rax_free(remove->supportVector);
151              rax_free(remove);
152              rax_free(tmp);
153            }
154          entryCount++;
155        }
156 
157       if(!start)
158         {
159           assert(!lastValid);
160           h->Items[k] = NULL;
161         }
162       else
163         {
164           h->Items[k] = start;
165         }
166     }
167 
168   assert(entryCount ==  h->entries);
169   h->entries-= removeCount;
170 }
171 
172 
173 
174 
175 
176 
177 
178 
179 
180 
181 
initBitVector(int mxtips,unsigned int * vectorLength)182 unsigned int **initBitVector(int mxtips, unsigned int *vectorLength)
183 {
184   unsigned int
185     **bitVectors = (unsigned int **)rax_malloc(sizeof(unsigned int*) * 2 * (size_t)mxtips);
186 
187   int
188     i;
189 
190   if(mxtips % PLL_MASK_LENGTH == 0)
191     *vectorLength = mxtips / PLL_MASK_LENGTH;
192   else
193     *vectorLength = 1 + (mxtips / PLL_MASK_LENGTH);
194 
195   for(i = 1; i <= mxtips; i++)
196     {
197       bitVectors[i] = (unsigned int *)rax_calloc((size_t)(*vectorLength), sizeof(unsigned int));
198       assert(bitVectors[i]);
199       bitVectors[i][(i - 1) / PLL_MASK_LENGTH] |= mask32[(i - 1) % PLL_MASK_LENGTH];
200     }
201 
202   for(i = mxtips + 1; i < 2 * mxtips; i++)
203     {
204       bitVectors[i] = (unsigned int *)rax_malloc(sizeof(unsigned int) * (size_t)(*vectorLength));
205       assert(bitVectors[i]);
206     }
207 
208   return bitVectors;
209 }
210 
freeBitVectors(unsigned int ** v,int n)211 void freeBitVectors(unsigned int **v, int n)
212 {
213   int i;
214 
215   for(i = 1; i < n; i++)
216     rax_free(v[i]);
217 }
218 
219 
newviewBipartitions(unsigned int ** bitVectors,nodeptr p,int numsp,unsigned int vectorLength,int processID)220 static void newviewBipartitions(unsigned int **bitVectors,
221                                 nodeptr p,
222                                 int numsp,
223                                 unsigned int vectorLength,
224                                 int processID)
225 {
226 
227   if(isTip(p->number, numsp))
228     return;
229   {
230     nodeptr
231       q = p->next->back,
232       r = p->next->next->back;
233 
234 
235 
236     unsigned int
237       *vector = bitVectors[p->number],
238       *left  = bitVectors[q->number],
239       *right = bitVectors[r->number];
240     unsigned
241       int i;
242 
243     assert(processID == 0);
244 
245 
246     while(!p->xBips)
247       {
248         if(!p->xBips)
249           getxnodeBips(p);
250       }
251 
252     p->hash = q->hash ^ r->hash;
253 
254     if(isTip(q->number, numsp) && isTip(r->number, numsp))
255       {
256         for(i = 0; i < vectorLength; i++)
257           vector[i] = left[i] | right[i];
258       }
259     else
260       {
261         if(isTip(q->number, numsp) || isTip(r->number, numsp))
262           {
263             if(isTip(r->number, numsp))
264               {
265                 nodeptr tmp = r;
266                 r = q;
267                 q = tmp;
268               }
269 
270             while(!r->xBips)
271               {
272                 if(!r->xBips)
273                   newviewBipartitions(bitVectors, r, numsp, vectorLength, processID);
274               }
275 
276             for(i = 0; i < vectorLength; i++)
277               vector[i] = left[i] | right[i];
278           }
279         else
280           {
281             while((!r->xBips) || (!q->xBips))
282               {
283                 if(!q->xBips)
284                   newviewBipartitions(bitVectors, q, numsp, vectorLength, processID);
285                 if(!r->xBips)
286                   newviewBipartitions(bitVectors, r, numsp, vectorLength, processID);
287               }
288 
289             for(i = 0; i < vectorLength; i++)
290               vector[i] = left[i] | right[i];
291           }
292 
293       }
294   }
295 }
296 
297 
298 
299 
insertHashRF(unsigned int * bitVector,pllHashTable * h,unsigned int vectorLength,int treeNumber,int treeVectorLength,hashNumberType position,int support,pllBoolean computeWRF)300 static void insertHashRF(unsigned int *bitVector,
301                          pllHashTable *h,
302                          unsigned int vectorLength,
303                          int treeNumber,
304                          int treeVectorLength,
305                          hashNumberType position,
306                          int support,
307                          pllBoolean computeWRF)
308 {
309   pllBipartitionEntry * e;
310   pllHashItem * hitem;
311 
312   if(h->Items[position] != NULL)
313     {
314       for (hitem = h->Items[position]; hitem; hitem = hitem->next)
315         {
316           e = (pllBipartitionEntry *)(hitem->data);
317 
318           if (!memcmp(bitVector, e->bitVector, vectorLength * sizeof(unsigned int)))
319             {
320               e->treeVector[treeNumber / PLL_MASK_LENGTH] |= mask32[treeNumber % PLL_MASK_LENGTH];
321               if(computeWRF)
322                 {
323                   e->supportVector[treeNumber] = support;
324                   assert(0 <= treeNumber && treeNumber < treeVectorLength * PLL_MASK_LENGTH);
325                 }
326               return;
327             }
328         }
329     }
330   e = initEntry();
331 
332   rax_posix_memalign ((void **)&(e->bitVector), PLL_BYTE_ALIGNMENT, (size_t)vectorLength * sizeof(unsigned int));
333   memset(e->bitVector, 0, vectorLength * sizeof(unsigned int));
334 
335   e->treeVector = (unsigned int*)rax_calloc((size_t)treeVectorLength, sizeof(unsigned int));
336   if(computeWRF)
337     e->supportVector = (int*)rax_calloc((size_t)treeVectorLength * PLL_MASK_LENGTH, sizeof(int));
338 
339   e->treeVector[treeNumber / PLL_MASK_LENGTH] |= mask32[treeNumber % PLL_MASK_LENGTH];
340   if(computeWRF)
341     {
342       e->supportVector[treeNumber] = support;
343 
344       assert(0 <= treeNumber && treeNumber < treeVectorLength * PLL_MASK_LENGTH);
345     }
346 
347   memcpy(e->bitVector, bitVector, sizeof(unsigned int) * vectorLength);
348 
349   pllHashAdd (h, position, NULL, (void *)e);
350 }
351 
352 
353 
bitVectorInitravSpecial(unsigned int ** bitVectors,nodeptr p,int numsp,unsigned int vectorLength,pllHashTable * h,int treeNumber,int function,branchInfo * bInf,int * countBranches,int treeVectorLength,pllBoolean traverseOnly,pllBoolean computeWRF,int processID)354 void bitVectorInitravSpecial(unsigned int **bitVectors, nodeptr p, int numsp, unsigned int vectorLength, pllHashTable *h, int treeNumber, int function, branchInfo *bInf,
355                              int *countBranches, int treeVectorLength, pllBoolean traverseOnly, pllBoolean computeWRF, int processID)
356 {
357   if(isTip(p->number, numsp))
358     return;
359   else
360     {
361       nodeptr
362         q = p->next;
363 
364       do
365         {
366           bitVectorInitravSpecial(bitVectors, q->back, numsp, vectorLength, h, treeNumber, function, bInf, countBranches, treeVectorLength, traverseOnly, computeWRF, processID);
367           q = q->next;
368         }
369       while(q != p);
370 
371       newviewBipartitions(bitVectors, p, numsp, vectorLength, processID);
372 
373       assert(p->xBips);
374 
375       assert(!traverseOnly);
376 
377       if(!(isTip(p->back->number, numsp)))
378         {
379           unsigned int
380             *toInsert  = bitVectors[p->number];
381 
382           hashNumberType
383             position = p->hash % h->size;
384 
385           assert(!(toInsert[0] & 1));
386           assert(!computeWRF);
387 
388           switch(function)
389             {
390             case PLL_BIPARTITIONS_RF:
391               insertHashRF(toInsert, h, vectorLength, treeNumber, treeVectorLength, position, 0, computeWRF);
392               *countBranches =  *countBranches + 1;
393               break;
394             default:
395               assert(0);
396             }
397         }
398 
399     }
400 }
401 
convergenceCriterion(pllHashTable * h,int mxtips)402 double convergenceCriterion(pllHashTable *h, int mxtips)
403 {
404   int
405     rf = 0;
406 
407   unsigned int
408     k = 0,
409     entryCount = 0;
410 
411   double
412     rrf;
413 
414   pllHashItem * hitem;
415 
416   for(k = 0, entryCount = 0; k < h->size; k++)
417     {
418       for (hitem = h->Items[k]; hitem; hitem = hitem->next)
419        {
420          pllBipartitionEntry *e = hitem->data;
421          unsigned int *vector = e->treeVector;
422 
423          if(((vector[0] & 1) > 0) + ((vector[0] & 2) > 0) == 1)
424            rf++;
425 
426          entryCount++;
427          e = e->next;
428        }
429     }
430 
431   assert(entryCount == h->entries);
432   rrf = (double)rf/((double)(2 * (mxtips - 3)));
433   return rrf;
434 }
435