1 | /* RAxML-VI-HPC (version 2.2) a program for sequential and parallel estimation of phylogenetic trees |
---|
2 | * Copyright August 2006 by Alexandros Stamatakis |
---|
3 | * |
---|
4 | * Partially derived from |
---|
5 | * fastDNAml, a program for estimation of phylogenetic trees from sequences by Gary J. Olsen |
---|
6 | * |
---|
7 | * and |
---|
8 | * |
---|
9 | * Programs of the PHYLIP package by Joe Felsenstein. |
---|
10 | * |
---|
11 | * This program is free software; you may redistribute it and/or modify its |
---|
12 | * under the terms of the GNU General Public License as published by the Free |
---|
13 | * Software Foundation; either version 2 of the License, or (at your option) |
---|
14 | * any later version. |
---|
15 | * |
---|
16 | * This program is distributed in the hope that it will be useful, but |
---|
17 | * WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY |
---|
18 | * or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License |
---|
19 | * for more details. |
---|
20 | * |
---|
21 | * |
---|
22 | * For any other enquiries send an Email to Alexandros Stamatakis |
---|
23 | * Alexandros.Stamatakis@epfl.ch |
---|
24 | * |
---|
25 | * When publishing work that is based on the results from RAxML-VI-HPC please cite: |
---|
26 | * |
---|
27 | * Alexandros Stamatakis:"RAxML-VI-HPC: maximum likelihood-based phylogenetic analyses with thousands of taxa and mixed models". |
---|
28 | * Bioinformatics 2006; doi: 10.1093/bioinformatics/btl446 |
---|
29 | */ |
---|
30 | |
---|
31 | #ifndef WIN32 |
---|
32 | #include <sys/times.h> |
---|
33 | #include <sys/types.h> |
---|
34 | #include <sys/time.h> |
---|
35 | #include <unistd.h> |
---|
36 | #endif |
---|
37 | |
---|
38 | #include <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 | |
---|
46 | #include "axml.h" |
---|
47 | |
---|
48 | extern char workdir[1024]; |
---|
49 | extern char run_id[128]; |
---|
50 | |
---|
51 | extern const char binaryStateNames[2]; |
---|
52 | extern const char dnaStateNames[4]; |
---|
53 | extern const char protStateNames[20]; |
---|
54 | extern const char genericStateNames[32]; |
---|
55 | |
---|
56 | static char getStateCharacter(int dataType, int state) |
---|
57 | { |
---|
58 | char |
---|
59 | result; |
---|
60 | |
---|
61 | assert(MIN_MODEL < dataType && dataType < MAX_MODEL); |
---|
62 | |
---|
63 | switch(dataType) |
---|
64 | { |
---|
65 | case BINARY_DATA: |
---|
66 | result = binaryStateNames[state]; |
---|
67 | break; |
---|
68 | case DNA_DATA: |
---|
69 | result = dnaStateNames[state]; |
---|
70 | break; |
---|
71 | case AA_DATA: |
---|
72 | result = protStateNames[state]; |
---|
73 | break; |
---|
74 | case GENERIC_32: |
---|
75 | result = genericStateNames[state]; |
---|
76 | break; |
---|
77 | default: |
---|
78 | assert(0); |
---|
79 | } |
---|
80 | |
---|
81 | return result; |
---|
82 | } |
---|
83 | |
---|
84 | static void makeP_Flex_Ancestral(double *rptr, double *EI, double *EIGN, int numberOfCategories, double *left, const int numStates) |
---|
85 | { |
---|
86 | int |
---|
87 | i, |
---|
88 | j, |
---|
89 | k; |
---|
90 | |
---|
91 | const int |
---|
92 | rates = numStates - 1, |
---|
93 | statesSquare = numStates * numStates; |
---|
94 | |
---|
95 | double |
---|
96 | z1 = 0.0, |
---|
97 | lz1[64], |
---|
98 | d1[64]; |
---|
99 | |
---|
100 | assert(numStates <= 64); |
---|
101 | |
---|
102 | for(i = 0; i < rates; i++) |
---|
103 | lz1[i] = EIGN[i] * z1; |
---|
104 | |
---|
105 | |
---|
106 | for(i = 0; i < numberOfCategories; i++) |
---|
107 | { |
---|
108 | for(j = 0; j < rates; j++) |
---|
109 | d1[j] = EXP (rptr[i] * lz1[j]); |
---|
110 | |
---|
111 | for(j = 0; j < numStates; j++) |
---|
112 | { |
---|
113 | left[statesSquare * i + numStates * j] = 1.0; |
---|
114 | |
---|
115 | for(k = 1; k < numStates; k++) |
---|
116 | left[statesSquare * i + numStates * j + k] = d1[k-1] * EI[rates * j + (k-1)]; |
---|
117 | } |
---|
118 | } |
---|
119 | } |
---|
120 | |
---|
121 | static void ancestralCat(double *v, double *sumBuffer, double *diagptable, int i, int numStates) |
---|
122 | { |
---|
123 | int |
---|
124 | l, |
---|
125 | j; |
---|
126 | |
---|
127 | double |
---|
128 | *ancestral = &sumBuffer[numStates * i], |
---|
129 | sum = 0.0, |
---|
130 | *term = (double*)rax_malloc(sizeof(double) * numStates); |
---|
131 | |
---|
132 | for(l = 0; l < numStates; l++) |
---|
133 | { |
---|
134 | double |
---|
135 | ump_x1 = 0.0; |
---|
136 | |
---|
137 | for(j = 0; j < numStates; j++) |
---|
138 | ump_x1 += v[j] * diagptable[l * numStates + j]; |
---|
139 | |
---|
140 | sum += ump_x1; |
---|
141 | term[l] = ump_x1; |
---|
142 | |
---|
143 | } |
---|
144 | |
---|
145 | for(l = 0; l < numStates; l++) |
---|
146 | ancestral[l] = term[l] / sum; |
---|
147 | |
---|
148 | rax_free(term); |
---|
149 | } |
---|
150 | |
---|
151 | static void newviewFlexCat_Ancestral(int tipCase, double *extEV, |
---|
152 | int *cptr, |
---|
153 | double *x1, double *x2, double *tipVector, |
---|
154 | unsigned char *tipX1, unsigned char *tipX2, |
---|
155 | int n, double *left, double *right, const int numStates, double *diagptable, double *sumBuffer) |
---|
156 | { |
---|
157 | double |
---|
158 | *x3 = (double *)rax_malloc(sizeof(double) * numStates), |
---|
159 | *le, *ri, *v, *vl, *vr, |
---|
160 | ump_x1, ump_x2, x1px2; |
---|
161 | int |
---|
162 | i, l, j, scale; |
---|
163 | |
---|
164 | const int |
---|
165 | statesSquare = numStates * numStates; |
---|
166 | |
---|
167 | switch(tipCase) |
---|
168 | { |
---|
169 | case TIP_TIP: |
---|
170 | { |
---|
171 | for (i = 0; i < n; i++) |
---|
172 | { |
---|
173 | le = &left[cptr[i] * statesSquare]; |
---|
174 | ri = &right[cptr[i] * statesSquare]; |
---|
175 | |
---|
176 | vl = &(tipVector[numStates * tipX1[i]]); |
---|
177 | vr = &(tipVector[numStates * tipX2[i]]); |
---|
178 | |
---|
179 | v = x3; |
---|
180 | |
---|
181 | for(l = 0; l < numStates; l++) |
---|
182 | v[l] = 0.0; |
---|
183 | |
---|
184 | for(l = 0; l < numStates; l++) |
---|
185 | { |
---|
186 | ump_x1 = 0.0; |
---|
187 | ump_x2 = 0.0; |
---|
188 | |
---|
189 | for(j = 0; j < numStates; j++) |
---|
190 | { |
---|
191 | ump_x1 += vl[j] * le[l * numStates + j]; |
---|
192 | ump_x2 += vr[j] * ri[l * numStates + j]; |
---|
193 | } |
---|
194 | |
---|
195 | x1px2 = ump_x1 * ump_x2; |
---|
196 | |
---|
197 | for(j = 0; j < numStates; j++) |
---|
198 | v[j] += x1px2 * extEV[l * numStates + j]; |
---|
199 | } |
---|
200 | |
---|
201 | ancestralCat(v, sumBuffer, &diagptable[cptr[i] * statesSquare], i, numStates); |
---|
202 | } |
---|
203 | } |
---|
204 | break; |
---|
205 | case TIP_INNER: |
---|
206 | { |
---|
207 | for (i = 0; i < n; i++) |
---|
208 | { |
---|
209 | le = &left[cptr[i] * statesSquare]; |
---|
210 | ri = &right[cptr[i] * statesSquare]; |
---|
211 | |
---|
212 | vl = &(tipVector[numStates * tipX1[i]]); |
---|
213 | vr = &x2[numStates * i]; |
---|
214 | v = x3; |
---|
215 | |
---|
216 | for(l = 0; l < numStates; l++) |
---|
217 | v[l] = 0.0; |
---|
218 | |
---|
219 | for(l = 0; l < numStates; l++) |
---|
220 | { |
---|
221 | ump_x1 = 0.0; |
---|
222 | ump_x2 = 0.0; |
---|
223 | |
---|
224 | for(j = 0; j < numStates; j++) |
---|
225 | { |
---|
226 | ump_x1 += vl[j] * le[l * numStates + j]; |
---|
227 | ump_x2 += vr[j] * ri[l * numStates + j]; |
---|
228 | } |
---|
229 | |
---|
230 | x1px2 = ump_x1 * ump_x2; |
---|
231 | |
---|
232 | for(j = 0; j < numStates; j++) |
---|
233 | v[j] += x1px2 * extEV[l * numStates + j]; |
---|
234 | } |
---|
235 | |
---|
236 | scale = 1; |
---|
237 | for(l = 0; scale && (l < numStates); l++) |
---|
238 | scale = ((v[l] < minlikelihood) && (v[l] > minusminlikelihood)); |
---|
239 | |
---|
240 | if(scale) |
---|
241 | { |
---|
242 | for(l = 0; l < numStates; l++) |
---|
243 | v[l] *= twotothe256; |
---|
244 | } |
---|
245 | |
---|
246 | ancestralCat(v, sumBuffer, &diagptable[cptr[i] * statesSquare], i, numStates); |
---|
247 | } |
---|
248 | } |
---|
249 | break; |
---|
250 | case INNER_INNER: |
---|
251 | for(i = 0; i < n; i++) |
---|
252 | { |
---|
253 | le = &left[cptr[i] * statesSquare]; |
---|
254 | ri = &right[cptr[i] * statesSquare]; |
---|
255 | |
---|
256 | vl = &x1[numStates * i]; |
---|
257 | vr = &x2[numStates * i]; |
---|
258 | v = x3; |
---|
259 | |
---|
260 | for(l = 0; l < numStates; l++) |
---|
261 | v[l] = 0.0; |
---|
262 | |
---|
263 | for(l = 0; l < numStates; l++) |
---|
264 | { |
---|
265 | ump_x1 = 0.0; |
---|
266 | ump_x2 = 0.0; |
---|
267 | |
---|
268 | for(j = 0; j < numStates; j++) |
---|
269 | { |
---|
270 | ump_x1 += vl[j] * le[l * numStates + j]; |
---|
271 | ump_x2 += vr[j] * ri[l * numStates + j]; |
---|
272 | } |
---|
273 | |
---|
274 | x1px2 = ump_x1 * ump_x2; |
---|
275 | |
---|
276 | for(j = 0; j < numStates; j++) |
---|
277 | v[j] += x1px2 * extEV[l * numStates + j]; |
---|
278 | } |
---|
279 | |
---|
280 | scale = 1; |
---|
281 | for(l = 0; scale && (l < numStates); l++) |
---|
282 | scale = ((v[l] < minlikelihood) && (v[l] > minusminlikelihood)); |
---|
283 | |
---|
284 | if(scale) |
---|
285 | { |
---|
286 | for(l = 0; l < numStates; l++) |
---|
287 | v[l] *= twotothe256; |
---|
288 | } |
---|
289 | |
---|
290 | ancestralCat(v, sumBuffer, &diagptable[cptr[i] * statesSquare], i, numStates); |
---|
291 | } |
---|
292 | break; |
---|
293 | default: |
---|
294 | assert(0); |
---|
295 | } |
---|
296 | |
---|
297 | rax_free(x3); |
---|
298 | |
---|
299 | |
---|
300 | |
---|
301 | } |
---|
302 | |
---|
303 | |
---|
304 | |
---|
305 | static void ancestralGamma(double *_v, double *sumBuffer, double *diagptable, int i, int numStates, int gammaStates) |
---|
306 | { |
---|
307 | int |
---|
308 | statesSquare = numStates * numStates, |
---|
309 | k, |
---|
310 | j, |
---|
311 | l; |
---|
312 | |
---|
313 | double |
---|
314 | *v, |
---|
315 | *ancestral = &sumBuffer[gammaStates * i], |
---|
316 | sum = 0.0, |
---|
317 | *term = (double*)rax_malloc(sizeof(double) * numStates); |
---|
318 | |
---|
319 | for(l = 0; l < numStates; l++) |
---|
320 | term[l] = 0.0; |
---|
321 | |
---|
322 | for(k = 0; k < 4; k++) |
---|
323 | { |
---|
324 | v = &(_v[numStates * k]); |
---|
325 | |
---|
326 | for(l = 0; l < numStates; l++) |
---|
327 | { |
---|
328 | double |
---|
329 | al = 0.0; |
---|
330 | |
---|
331 | for(j = 0; j < numStates; j++) |
---|
332 | al += v[j] * diagptable[k * statesSquare + l * numStates + j]; |
---|
333 | |
---|
334 | term[l] += al; |
---|
335 | sum += al; |
---|
336 | } |
---|
337 | } |
---|
338 | |
---|
339 | for(l = 0; l < numStates; l++) |
---|
340 | ancestral[l] = term[l] / sum; |
---|
341 | |
---|
342 | rax_free(term); |
---|
343 | } |
---|
344 | |
---|
345 | |
---|
346 | static void newviewFlexGamma_Ancestral(int tipCase, |
---|
347 | double *x1, double *x2, double *extEV, double *tipVector, |
---|
348 | unsigned char *tipX1, unsigned char *tipX2, |
---|
349 | int n, double *left, double *right, |
---|
350 | const int numStates, double *diagptable, double *sumBuffer) |
---|
351 | { |
---|
352 | double |
---|
353 | *v, |
---|
354 | *x3 = (double*)rax_malloc(sizeof(double) * 4 * numStates); |
---|
355 | double x1px2; |
---|
356 | int i, j, l, k, scale; |
---|
357 | double *vl, *vr, al, ar; |
---|
358 | |
---|
359 | |
---|
360 | |
---|
361 | const int |
---|
362 | statesSquare = numStates * numStates, |
---|
363 | gammaStates = 4 * numStates; |
---|
364 | |
---|
365 | switch(tipCase) |
---|
366 | { |
---|
367 | case TIP_TIP: |
---|
368 | { |
---|
369 | for(i = 0; i < n; i++) |
---|
370 | { |
---|
371 | for(k = 0; k < 4; k++) |
---|
372 | { |
---|
373 | vl = &(tipVector[numStates * tipX1[i]]); |
---|
374 | vr = &(tipVector[numStates * tipX2[i]]); |
---|
375 | v = &(x3[numStates * k]); |
---|
376 | |
---|
377 | for(l = 0; l < numStates; l++) |
---|
378 | v[l] = 0; |
---|
379 | |
---|
380 | for(l = 0; l < numStates; l++) |
---|
381 | { |
---|
382 | al = 0.0; |
---|
383 | ar = 0.0; |
---|
384 | for(j = 0; j < numStates; j++) |
---|
385 | { |
---|
386 | al += vl[j] * left[k * statesSquare + l * numStates + j]; |
---|
387 | ar += vr[j] * right[k * statesSquare + l * numStates + j]; |
---|
388 | } |
---|
389 | |
---|
390 | x1px2 = al * ar; |
---|
391 | for(j = 0; j < numStates; j++) |
---|
392 | v[j] += x1px2 * extEV[numStates * l + j]; |
---|
393 | } |
---|
394 | } |
---|
395 | |
---|
396 | ancestralGamma(x3, sumBuffer, diagptable, i, numStates, gammaStates); |
---|
397 | } |
---|
398 | } |
---|
399 | break; |
---|
400 | case TIP_INNER: |
---|
401 | { |
---|
402 | for (i = 0; i < n; i++) |
---|
403 | { |
---|
404 | for(k = 0; k < 4; k++) |
---|
405 | { |
---|
406 | vl = &(tipVector[numStates * tipX1[i]]); |
---|
407 | vr = &(x2[gammaStates * i + numStates * k]); |
---|
408 | v = &(x3[numStates * k]); |
---|
409 | |
---|
410 | for(l = 0; l < numStates; l++) |
---|
411 | v[l] = 0; |
---|
412 | |
---|
413 | for(l = 0; l < numStates; l++) |
---|
414 | { |
---|
415 | al = 0.0; |
---|
416 | ar = 0.0; |
---|
417 | for(j = 0; j < numStates; j++) |
---|
418 | { |
---|
419 | al += vl[j] * left[k * statesSquare + l * numStates + j]; |
---|
420 | ar += vr[j] * right[k * statesSquare + l * numStates + j]; |
---|
421 | } |
---|
422 | |
---|
423 | x1px2 = al * ar; |
---|
424 | for(j = 0; j < numStates; j++) |
---|
425 | v[j] += x1px2 * extEV[numStates * l + j]; |
---|
426 | } |
---|
427 | } |
---|
428 | |
---|
429 | v = x3; |
---|
430 | scale = 1; |
---|
431 | for(l = 0; scale && (l < gammaStates); l++) |
---|
432 | scale = (ABS(v[l]) < minlikelihood); |
---|
433 | |
---|
434 | if(scale) |
---|
435 | { |
---|
436 | for(l = 0; l < gammaStates; l++) |
---|
437 | v[l] *= twotothe256; |
---|
438 | } |
---|
439 | |
---|
440 | ancestralGamma(x3, sumBuffer, diagptable, i, numStates, gammaStates); |
---|
441 | } |
---|
442 | } |
---|
443 | break; |
---|
444 | case INNER_INNER: |
---|
445 | for (i = 0; i < n; i++) |
---|
446 | { |
---|
447 | for(k = 0; k < 4; k++) |
---|
448 | { |
---|
449 | vl = &(x1[gammaStates * i + numStates * k]); |
---|
450 | vr = &(x2[gammaStates * i + numStates * k]); |
---|
451 | v = &(x3[numStates * k]); |
---|
452 | |
---|
453 | for(l = 0; l < numStates; l++) |
---|
454 | v[l] = 0; |
---|
455 | |
---|
456 | for(l = 0; l < numStates; l++) |
---|
457 | { |
---|
458 | al = 0.0; |
---|
459 | ar = 0.0; |
---|
460 | for(j = 0; j < numStates; j++) |
---|
461 | { |
---|
462 | al += vl[j] * left[k * statesSquare + l * numStates + j]; |
---|
463 | ar += vr[j] * right[k * statesSquare + l * numStates + j]; |
---|
464 | } |
---|
465 | |
---|
466 | x1px2 = al * ar; |
---|
467 | for(j = 0; j < numStates; j++) |
---|
468 | v[j] += x1px2 * extEV[numStates * l + j]; |
---|
469 | } |
---|
470 | } |
---|
471 | |
---|
472 | v = x3; |
---|
473 | scale = 1; |
---|
474 | for(l = 0; scale && (l < gammaStates); l++) |
---|
475 | scale = ((ABS(v[l]) < minlikelihood)); |
---|
476 | |
---|
477 | if (scale) |
---|
478 | { |
---|
479 | for(l = 0; l < gammaStates; l++) |
---|
480 | v[l] *= twotothe256; |
---|
481 | } |
---|
482 | |
---|
483 | ancestralGamma(x3, sumBuffer, diagptable, i, numStates, gammaStates); |
---|
484 | } |
---|
485 | break; |
---|
486 | default: |
---|
487 | assert(0); |
---|
488 | } |
---|
489 | |
---|
490 | |
---|
491 | |
---|
492 | rax_free(x3); |
---|
493 | |
---|
494 | } |
---|
495 | void newviewIterativeAncestral(tree *tr) |
---|
496 | { |
---|
497 | traversalInfo |
---|
498 | *ti = tr->td[0].ti; |
---|
499 | |
---|
500 | int |
---|
501 | i, |
---|
502 | model; |
---|
503 | |
---|
504 | assert(!tr->saveMemory); |
---|
505 | |
---|
506 | for(i = 1; i < tr->td[0].count; i++) |
---|
507 | { |
---|
508 | traversalInfo |
---|
509 | *tInfo = &ti[i]; |
---|
510 | |
---|
511 | for(model = 0; model < tr->NumberOfModels; model++) |
---|
512 | { |
---|
513 | double |
---|
514 | *x1_start = (double*)NULL, |
---|
515 | *x2_start = (double*)NULL, |
---|
516 | *left = (double*)NULL, |
---|
517 | *right = (double*)NULL, |
---|
518 | qz, |
---|
519 | rz; |
---|
520 | |
---|
521 | unsigned char |
---|
522 | *tipX1 = (unsigned char *)NULL, |
---|
523 | *tipX2 = (unsigned char *)NULL; |
---|
524 | |
---|
525 | size_t |
---|
526 | states = (size_t)tr->partitionData[model].states, |
---|
527 | width = tr->partitionData[model].width; |
---|
528 | |
---|
529 | |
---|
530 | |
---|
531 | switch(tInfo->tipCase) |
---|
532 | { |
---|
533 | case TIP_TIP: |
---|
534 | tipX1 = tr->partitionData[model].yVector[tInfo->qNumber]; |
---|
535 | tipX2 = tr->partitionData[model].yVector[tInfo->rNumber]; |
---|
536 | break; |
---|
537 | case TIP_INNER: |
---|
538 | tipX1 = tr->partitionData[model].yVector[tInfo->qNumber]; |
---|
539 | x2_start = tr->partitionData[model].xVector[tInfo->rNumber - tr->mxtips - 1]; |
---|
540 | break; |
---|
541 | case INNER_INNER: |
---|
542 | x1_start = tr->partitionData[model].xVector[tInfo->qNumber - tr->mxtips - 1]; |
---|
543 | x2_start = tr->partitionData[model].xVector[tInfo->rNumber - tr->mxtips - 1]; |
---|
544 | break; |
---|
545 | default: |
---|
546 | assert(0); |
---|
547 | } |
---|
548 | |
---|
549 | left = tr->partitionData[model].left; |
---|
550 | right = tr->partitionData[model].right; |
---|
551 | |
---|
552 | if(tr->multiBranch) |
---|
553 | { |
---|
554 | qz = tInfo->qz[model]; |
---|
555 | rz = tInfo->rz[model]; |
---|
556 | } |
---|
557 | else |
---|
558 | { |
---|
559 | qz = tInfo->qz[0]; |
---|
560 | rz = tInfo->rz[0]; |
---|
561 | } |
---|
562 | |
---|
563 | switch(tr->rateHetModel) |
---|
564 | { |
---|
565 | case CAT: |
---|
566 | { |
---|
567 | double |
---|
568 | *diagptable = (double*)rax_malloc(tr->partitionData[model].numberOfCategories * states * states * sizeof(double)); |
---|
569 | |
---|
570 | makeP_Flex(qz, rz, tr->partitionData[model].perSiteRates, |
---|
571 | tr->partitionData[model].EI, |
---|
572 | tr->partitionData[model].EIGN, |
---|
573 | tr->partitionData[model].numberOfCategories, left, right, states); |
---|
574 | |
---|
575 | |
---|
576 | makeP_Flex_Ancestral(tr->partitionData[model].perSiteRates, |
---|
577 | tr->partitionData[model].EI, |
---|
578 | tr->partitionData[model].EIGN, |
---|
579 | tr->partitionData[model].numberOfCategories, diagptable, states); |
---|
580 | |
---|
581 | |
---|
582 | newviewFlexCat_Ancestral(tInfo->tipCase, tr->partitionData[model].EV, tr->partitionData[model].rateCategory, |
---|
583 | x1_start, x2_start, tr->partitionData[model].tipVector, |
---|
584 | tipX1, tipX2, width, left, right, states, diagptable, |
---|
585 | tr->partitionData[model].sumBuffer); |
---|
586 | |
---|
587 | rax_free(diagptable); |
---|
588 | } |
---|
589 | break; |
---|
590 | case GAMMA: |
---|
591 | case GAMMA_I: |
---|
592 | { |
---|
593 | double |
---|
594 | *diagptable = (double*)rax_malloc(4 * states * states * sizeof(double)); |
---|
595 | |
---|
596 | makeP_Flex(qz, rz, tr->partitionData[model].gammaRates, |
---|
597 | tr->partitionData[model].EI, |
---|
598 | tr->partitionData[model].EIGN, |
---|
599 | 4, left, right, states); |
---|
600 | |
---|
601 | makeP_Flex_Ancestral(tr->partitionData[model].gammaRates, |
---|
602 | tr->partitionData[model].EI, |
---|
603 | tr->partitionData[model].EIGN, |
---|
604 | 4, diagptable, states); |
---|
605 | |
---|
606 | |
---|
607 | newviewFlexGamma_Ancestral(tInfo->tipCase, |
---|
608 | x1_start, x2_start, |
---|
609 | tr->partitionData[model].EV, |
---|
610 | tr->partitionData[model].tipVector, |
---|
611 | tipX1, tipX2, |
---|
612 | width, left, right, states, diagptable, tr->partitionData[model].sumBuffer); |
---|
613 | |
---|
614 | rax_free(diagptable); |
---|
615 | } |
---|
616 | break; |
---|
617 | default: |
---|
618 | assert(0); |
---|
619 | } |
---|
620 | } |
---|
621 | } |
---|
622 | |
---|
623 | } |
---|
624 | static void traversalInfoAncestralRoot(nodeptr p, traversalInfo *ti, int *counter, int maxTips, int numBranches) |
---|
625 | { |
---|
626 | int |
---|
627 | i; |
---|
628 | |
---|
629 | nodeptr |
---|
630 | q = p, |
---|
631 | r = p->back; |
---|
632 | |
---|
633 | for(i = 0; i < numBranches; i++) |
---|
634 | { |
---|
635 | double |
---|
636 | z = sqrt(q->z[i]); |
---|
637 | if(z < zmin) |
---|
638 | z = zmin; |
---|
639 | if(z > zmax) |
---|
640 | z = zmax; |
---|
641 | |
---|
642 | z = log(z); |
---|
643 | |
---|
644 | ti[*counter].qz[i] = z; |
---|
645 | ti[*counter].rz[i] = z; |
---|
646 | } |
---|
647 | |
---|
648 | if(isTip(r->number, maxTips) && isTip(q->number, maxTips)) |
---|
649 | { |
---|
650 | ti[*counter].tipCase = TIP_TIP; |
---|
651 | ti[*counter].qNumber = q->number; |
---|
652 | ti[*counter].rNumber = r->number; |
---|
653 | |
---|
654 | *counter = *counter + 1; |
---|
655 | } |
---|
656 | else |
---|
657 | { |
---|
658 | if(isTip(r->number, maxTips) || isTip(q->number, maxTips)) |
---|
659 | { |
---|
660 | nodeptr tmp; |
---|
661 | |
---|
662 | if(isTip(r->number, maxTips)) |
---|
663 | { |
---|
664 | tmp = r; |
---|
665 | r = q; |
---|
666 | q = tmp; |
---|
667 | } |
---|
668 | |
---|
669 | ti[*counter].tipCase = TIP_INNER; |
---|
670 | ti[*counter].qNumber = q->number; |
---|
671 | ti[*counter].rNumber = r->number; |
---|
672 | |
---|
673 | *counter = *counter + 1; |
---|
674 | } |
---|
675 | else |
---|
676 | { |
---|
677 | ti[*counter].tipCase = INNER_INNER; |
---|
678 | ti[*counter].qNumber = q->number; |
---|
679 | ti[*counter].rNumber = r->number; |
---|
680 | |
---|
681 | *counter = *counter + 1; |
---|
682 | } |
---|
683 | } |
---|
684 | } |
---|
685 | |
---|
686 | |
---|
687 | void newviewGenericAncestral (tree *tr, nodeptr p, boolean atRoot) |
---|
688 | { |
---|
689 | if(atRoot) |
---|
690 | { |
---|
691 | tr->td[0].count = 1; |
---|
692 | traversalInfoAncestralRoot(p, &(tr->td[0].ti[0]), &(tr->td[0].count), tr->mxtips, tr->numBranches); |
---|
693 | |
---|
694 | if(tr->td[0].count > 1) |
---|
695 | { |
---|
696 | #ifdef _USE_PTHREADS |
---|
697 | masterBarrier(THREAD_NEWVIEW_ANCESTRAL, tr); |
---|
698 | #else |
---|
699 | newviewIterativeAncestral(tr); |
---|
700 | #endif |
---|
701 | } |
---|
702 | } |
---|
703 | else |
---|
704 | { |
---|
705 | if(isTip(p->number, tr->mxtips)) |
---|
706 | return; |
---|
707 | |
---|
708 | |
---|
709 | tr->td[0].count = 1; |
---|
710 | computeTraversalInfo(p, &(tr->td[0].ti[0]), &(tr->td[0].count), tr->mxtips, tr->numBranches); |
---|
711 | |
---|
712 | if(tr->td[0].count > 1) |
---|
713 | { |
---|
714 | #ifdef _USE_PTHREADS |
---|
715 | masterBarrier(THREAD_NEWVIEW_ANCESTRAL, tr); |
---|
716 | #else |
---|
717 | newviewIterativeAncestral(tr); |
---|
718 | #endif |
---|
719 | } |
---|
720 | } |
---|
721 | } |
---|
722 | |
---|
723 | typedef struct { |
---|
724 | double *probs; |
---|
725 | char c; |
---|
726 | int states; |
---|
727 | } ancestralState; |
---|
728 | |
---|
729 | |
---|
730 | static void computeAncestralRec(tree *tr, nodeptr p, int *counter, FILE *probsFile, FILE *statesFile, boolean atRoot) |
---|
731 | { |
---|
732 | #ifdef _USE_PTHREADS |
---|
733 | size_t |
---|
734 | accumulatedOffset = 0; |
---|
735 | #endif |
---|
736 | |
---|
737 | int |
---|
738 | model, |
---|
739 | globalIndex = 0; |
---|
740 | |
---|
741 | ancestralState |
---|
742 | *a = (ancestralState *)rax_malloc(sizeof(ancestralState) * tr->cdta->endsite), |
---|
743 | *unsortedA = (ancestralState *)rax_malloc(sizeof(ancestralState) * tr->rdta->sites); |
---|
744 | |
---|
745 | if(!atRoot) |
---|
746 | { |
---|
747 | if(isTip(p->number, tr->mxtips)) |
---|
748 | return; |
---|
749 | |
---|
750 | computeAncestralRec(tr, p->next->back, counter, probsFile, statesFile, atRoot); |
---|
751 | computeAncestralRec(tr, p->next->next->back, counter, probsFile, statesFile, atRoot); |
---|
752 | |
---|
753 | newviewGeneric(tr, p); |
---|
754 | } |
---|
755 | |
---|
756 | newviewGenericAncestral(tr, p, atRoot); |
---|
757 | |
---|
758 | #ifdef _USE_PTHREADS |
---|
759 | masterBarrier(THREAD_GATHER_ANCESTRAL, tr); |
---|
760 | #endif |
---|
761 | |
---|
762 | if(atRoot) |
---|
763 | { |
---|
764 | fprintf(probsFile, "ROOT\n"); |
---|
765 | fprintf(statesFile, "ROOT "); |
---|
766 | } |
---|
767 | else |
---|
768 | { |
---|
769 | fprintf(probsFile, "%d\n", p->number); |
---|
770 | fprintf(statesFile, "%d ", p->number); |
---|
771 | } |
---|
772 | |
---|
773 | for(model = 0; model < tr->NumberOfModels; model++) |
---|
774 | { |
---|
775 | int |
---|
776 | offset, |
---|
777 | i, |
---|
778 | width = tr->partitionData[model].upper - tr->partitionData[model].lower, |
---|
779 | states = tr->partitionData[model].states; |
---|
780 | #ifdef _USE_PTHREADS |
---|
781 | double |
---|
782 | *ancestral = &tr->ancestralStates[accumulatedOffset]; |
---|
783 | #else |
---|
784 | double |
---|
785 | *ancestral = tr->partitionData[model].sumBuffer; |
---|
786 | #endif |
---|
787 | |
---|
788 | if(tr->rateHetModel == CAT) |
---|
789 | offset = 1; |
---|
790 | else |
---|
791 | offset = 4; |
---|
792 | |
---|
793 | for(i = 0; i < width; i++, globalIndex++) |
---|
794 | { |
---|
795 | double |
---|
796 | equal = 1.0 / (double)states, |
---|
797 | max = -1.0; |
---|
798 | |
---|
799 | boolean |
---|
800 | approximatelyEqual = TRUE; |
---|
801 | |
---|
802 | int |
---|
803 | max_l = -1, |
---|
804 | l; |
---|
805 | |
---|
806 | char |
---|
807 | c; |
---|
808 | |
---|
809 | a[globalIndex].states = states; |
---|
810 | a[globalIndex].probs = (double *)rax_malloc(sizeof(double) * states); |
---|
811 | |
---|
812 | for(l = 0; l < states; l++) |
---|
813 | { |
---|
814 | double |
---|
815 | value = ancestral[offset * states * i + l]; |
---|
816 | |
---|
817 | if(value > max) |
---|
818 | { |
---|
819 | max = value; |
---|
820 | max_l = l; |
---|
821 | } |
---|
822 | |
---|
823 | approximatelyEqual = approximatelyEqual && (ABS(equal - value) < 0.000001); |
---|
824 | |
---|
825 | a[globalIndex].probs[l] = value; |
---|
826 | } |
---|
827 | |
---|
828 | |
---|
829 | if(approximatelyEqual) |
---|
830 | c = '?'; |
---|
831 | else |
---|
832 | c = getStateCharacter(tr->partitionData[model].dataType, max_l); |
---|
833 | |
---|
834 | a[globalIndex].c = c; |
---|
835 | } |
---|
836 | |
---|
837 | #ifdef _USE_PTHREADS |
---|
838 | accumulatedOffset += width * offset * states; |
---|
839 | #endif |
---|
840 | } |
---|
841 | |
---|
842 | { |
---|
843 | int |
---|
844 | j, |
---|
845 | k; |
---|
846 | |
---|
847 | for(j = 0; j < tr->cdta->endsite; j++) |
---|
848 | { |
---|
849 | for(k = 0; k < tr->rdta->sites; k++) |
---|
850 | if(j == tr->patternPosition[k]) |
---|
851 | { |
---|
852 | int |
---|
853 | sorted = j, |
---|
854 | unsorted = tr->columnPosition[k] - 1; |
---|
855 | |
---|
856 | unsortedA[unsorted].states = a[sorted].states; |
---|
857 | unsortedA[unsorted].c = a[sorted].c; |
---|
858 | unsortedA[unsorted].probs = (double*)rax_malloc(sizeof(double) * unsortedA[unsorted].states); |
---|
859 | memcpy(unsortedA[unsorted].probs, a[sorted].probs, sizeof(double) * a[sorted].states); |
---|
860 | } |
---|
861 | } |
---|
862 | |
---|
863 | for(k = 0; k < tr->rdta->sites; k++) |
---|
864 | { |
---|
865 | for(j = 0; j < unsortedA[k].states; j++) |
---|
866 | fprintf(probsFile, "%f ", unsortedA[k].probs[j]); |
---|
867 | fprintf(probsFile, "\n"); |
---|
868 | fprintf(statesFile, "%c", unsortedA[k].c); |
---|
869 | } |
---|
870 | fprintf(probsFile, "\n"); |
---|
871 | fprintf(statesFile, "\n"); |
---|
872 | } |
---|
873 | |
---|
874 | |
---|
875 | *counter = *counter + 1; |
---|
876 | |
---|
877 | { |
---|
878 | int j; |
---|
879 | |
---|
880 | for(j = 0; j < tr->rdta->sites; j++) |
---|
881 | rax_free(unsortedA[j].probs); |
---|
882 | for(j = 0; j < tr->cdta->endsite; j++) |
---|
883 | rax_free(a[j].probs); |
---|
884 | } |
---|
885 | |
---|
886 | rax_free(a); |
---|
887 | rax_free(unsortedA); |
---|
888 | } |
---|
889 | |
---|
890 | static char *ancestralTreeRec(char *treestr, tree *tr, nodeptr p) |
---|
891 | { |
---|
892 | if(isTip(p->number, tr->rdta->numsp)) |
---|
893 | { |
---|
894 | sprintf(treestr, "%s", tr->nameList[p->number]); |
---|
895 | |
---|
896 | while (*treestr) |
---|
897 | treestr++; |
---|
898 | } |
---|
899 | else |
---|
900 | { |
---|
901 | *treestr++ = '('; |
---|
902 | treestr = ancestralTreeRec(treestr, tr, p->next->back); |
---|
903 | *treestr++ = ','; |
---|
904 | treestr = ancestralTreeRec(treestr, tr, p->next->next->back); |
---|
905 | *treestr++ = ')'; |
---|
906 | sprintf(treestr, "%d", p->number); |
---|
907 | } |
---|
908 | |
---|
909 | while (*treestr) |
---|
910 | treestr++; |
---|
911 | |
---|
912 | return treestr; |
---|
913 | } |
---|
914 | |
---|
915 | |
---|
916 | static char *ancestralTree(char *treestr, tree *tr) |
---|
917 | { |
---|
918 | *treestr++ = '('; |
---|
919 | treestr = ancestralTreeRec(treestr, tr, tr->leftRootNode); |
---|
920 | *treestr++ = ','; |
---|
921 | treestr = ancestralTreeRec(treestr, tr, tr->rightRootNode); |
---|
922 | *treestr++ = ')'; |
---|
923 | sprintf(treestr, "ROOT"); |
---|
924 | while (*treestr) |
---|
925 | treestr++; |
---|
926 | |
---|
927 | *treestr++ = ';'; |
---|
928 | |
---|
929 | *treestr++ = '\0'; |
---|
930 | while (*treestr) treestr++; |
---|
931 | return treestr; |
---|
932 | } |
---|
933 | |
---|
934 | void computeAncestralStates(tree *tr, double referenceLikelihood) |
---|
935 | { |
---|
936 | int |
---|
937 | counter = 0; |
---|
938 | |
---|
939 | char |
---|
940 | treeFileName[2048], |
---|
941 | ancestralProbsFileName[2048], |
---|
942 | ancestralStatesFileName[2048]; |
---|
943 | |
---|
944 | FILE |
---|
945 | *treeFile, |
---|
946 | *probsFile, |
---|
947 | *statesFile; |
---|
948 | |
---|
949 | #ifdef _USE_PTHREADS |
---|
950 | tr->ancestralStates = (double*)rax_malloc(getContiguousVectorLength(tr) * sizeof(double)); |
---|
951 | #endif |
---|
952 | |
---|
953 | strcpy(ancestralProbsFileName, workdir); |
---|
954 | strcpy(ancestralStatesFileName, workdir); |
---|
955 | strcpy(treeFileName, workdir); |
---|
956 | |
---|
957 | strcat(ancestralProbsFileName, "RAxML_marginalAncestralProbabilities."); |
---|
958 | strcat(ancestralStatesFileName, "RAxML_marginalAncestralStates."); |
---|
959 | strcat(treeFileName, "RAxML_nodeLabelledRootedTree."); |
---|
960 | |
---|
961 | strcat(ancestralProbsFileName, run_id); |
---|
962 | strcat(ancestralStatesFileName, run_id); |
---|
963 | strcat(treeFileName, run_id); |
---|
964 | |
---|
965 | probsFile = myfopen(ancestralProbsFileName, "w"); |
---|
966 | statesFile = myfopen(ancestralStatesFileName, "w"); |
---|
967 | treeFile = myfopen(treeFileName, "w"); |
---|
968 | |
---|
969 | assert(tr->leftRootNode == tr->rightRootNode->back); |
---|
970 | |
---|
971 | computeAncestralRec(tr, tr->leftRootNode, &counter, probsFile, statesFile, FALSE); |
---|
972 | computeAncestralRec(tr, tr->rightRootNode, &counter, probsFile, statesFile, FALSE); |
---|
973 | computeAncestralRec(tr, tr->rightRootNode, &counter, probsFile, statesFile, TRUE); |
---|
974 | |
---|
975 | evaluateGeneric(tr, tr->rightRootNode); |
---|
976 | |
---|
977 | if(fabs(tr->likelihood - referenceLikelihood) > 0.5) |
---|
978 | { |
---|
979 | printf("Something suspiciuous is going on with the marginal ancestral probability computations\n"); |
---|
980 | assert(0); |
---|
981 | } |
---|
982 | |
---|
983 | assert(counter == tr->mxtips - 1); |
---|
984 | |
---|
985 | ancestralTree(tr->tree_string, tr); |
---|
986 | |
---|
987 | fprintf(treeFile, "%s\n", tr->tree_string); |
---|
988 | |
---|
989 | fclose(probsFile); |
---|
990 | fclose(statesFile); |
---|
991 | fclose(treeFile); |
---|
992 | |
---|
993 | printBothOpen("Marginal Ancestral Probabilities written to file:\n%s\n\n", ancestralProbsFileName); |
---|
994 | printBothOpen("Ancestral Sequences based on Marginal Ancestral Probabilities written to file:\n%s\n\n", ancestralStatesFileName); |
---|
995 | printBothOpen("Node-laballed ROOTED tree written to file:\n%s\n", treeFileName); |
---|
996 | } |
---|
997 | |
---|
998 | |
---|