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