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 <unistd.h> |
---|
33 | #endif |
---|
34 | |
---|
35 | #include <math.h> |
---|
36 | #include <time.h> |
---|
37 | #include <stdlib.h> |
---|
38 | #include <stdio.h> |
---|
39 | #include <ctype.h> |
---|
40 | #include <string.h> |
---|
41 | #include "axml.h" |
---|
42 | |
---|
43 | #ifdef __SIM_SSE3 |
---|
44 | #include <xmmintrin.h> |
---|
45 | #include <pmmintrin.h> |
---|
46 | #endif |
---|
47 | |
---|
48 | |
---|
49 | /********************** GTRCAT ***************************************/ |
---|
50 | |
---|
51 | |
---|
52 | #ifdef __SIM_SSE3 |
---|
53 | |
---|
54 | static inline void computeVectorGTRCATPROT(double *lVector, int *eVector, double ki, int i, double qz, double rz, |
---|
55 | traversalInfo *ti, double *EIGN, double *EI, double *EV, double *tipVector, |
---|
56 | unsigned char **yVector, int mxtips) |
---|
57 | { |
---|
58 | double *x1, *x2, *x3; |
---|
59 | int |
---|
60 | pNumber = ti->pNumber, |
---|
61 | rNumber = ti->rNumber, |
---|
62 | qNumber = ti->qNumber; |
---|
63 | |
---|
64 | x3 = &(lVector[20 * (pNumber - mxtips)]); |
---|
65 | |
---|
66 | switch(ti->tipCase) |
---|
67 | { |
---|
68 | case TIP_TIP: |
---|
69 | x1 = &(tipVector[20 * yVector[qNumber][i]]); |
---|
70 | x2 = &(tipVector[20 * yVector[rNumber][i]]); |
---|
71 | break; |
---|
72 | case TIP_INNER: |
---|
73 | x1 = &(tipVector[20 * yVector[qNumber][i]]); |
---|
74 | x2 = &( lVector[20 * (rNumber - mxtips)]); |
---|
75 | break; |
---|
76 | case INNER_INNER: |
---|
77 | x1 = &(lVector[20 * (qNumber - mxtips)]); |
---|
78 | x2 = &(lVector[20 * (rNumber - mxtips)]); |
---|
79 | break; |
---|
80 | default: |
---|
81 | assert(0); |
---|
82 | } |
---|
83 | |
---|
84 | { |
---|
85 | double |
---|
86 | e1[20] __attribute__ ((aligned (BYTE_ALIGNMENT))), |
---|
87 | e2[20] __attribute__ ((aligned (BYTE_ALIGNMENT))), |
---|
88 | d1[20] __attribute__ ((aligned (BYTE_ALIGNMENT))), |
---|
89 | d2[20] __attribute__ ((aligned (BYTE_ALIGNMENT))), |
---|
90 | lz1, lz2; |
---|
91 | int l, k, scale; |
---|
92 | |
---|
93 | lz1 = qz * ki; |
---|
94 | lz2 = rz * ki; |
---|
95 | |
---|
96 | e1[0] = 1.0; |
---|
97 | e2[0] = 1.0; |
---|
98 | |
---|
99 | for(l = 1; l < 20; l++) |
---|
100 | { |
---|
101 | e1[l] = EXP(EIGN[l - 1] * lz1); |
---|
102 | e2[l] = EXP(EIGN[l - 1] * lz2); |
---|
103 | } |
---|
104 | |
---|
105 | for(l = 0; l < 20; l+=2) |
---|
106 | { |
---|
107 | __m128d d1v = _mm_mul_pd(_mm_load_pd(&x1[l]), _mm_load_pd(&e1[l])); |
---|
108 | __m128d d2v = _mm_mul_pd(_mm_load_pd(&x2[l]), _mm_load_pd(&e2[l])); |
---|
109 | |
---|
110 | _mm_store_pd(&d1[l], d1v); |
---|
111 | _mm_store_pd(&d2[l], d2v); |
---|
112 | } |
---|
113 | |
---|
114 | __m128d zero = _mm_setzero_pd(); |
---|
115 | |
---|
116 | for(l = 0; l < 20; l+=2) |
---|
117 | _mm_store_pd(&x3[l], zero); |
---|
118 | |
---|
119 | for(l = 0; l < 20; l++) |
---|
120 | { |
---|
121 | double *ev = &EV[l * 20]; |
---|
122 | __m128d ump_x1v = _mm_setzero_pd(); |
---|
123 | __m128d ump_x2v = _mm_setzero_pd(); |
---|
124 | __m128d x1px2v; |
---|
125 | |
---|
126 | for(k = 0; k < 20; k+=2) |
---|
127 | { |
---|
128 | __m128d eiv = _mm_load_pd(&EI[20 * l + k]); |
---|
129 | __m128d d1v = _mm_load_pd(&d1[k]); |
---|
130 | __m128d d2v = _mm_load_pd(&d2[k]); |
---|
131 | ump_x1v = _mm_add_pd(ump_x1v, _mm_mul_pd(d1v, eiv)); |
---|
132 | ump_x2v = _mm_add_pd(ump_x2v, _mm_mul_pd(d2v, eiv)); |
---|
133 | } |
---|
134 | |
---|
135 | ump_x1v = _mm_hadd_pd(ump_x1v, ump_x1v); |
---|
136 | ump_x2v = _mm_hadd_pd(ump_x2v, ump_x2v); |
---|
137 | |
---|
138 | x1px2v = _mm_mul_pd(ump_x1v, ump_x2v); |
---|
139 | |
---|
140 | for(k = 0; k < 20; k+=2) |
---|
141 | { |
---|
142 | __m128d ex3v = _mm_load_pd(&x3[k]); |
---|
143 | __m128d EVV = _mm_load_pd(&ev[k]); |
---|
144 | ex3v = _mm_add_pd(ex3v, _mm_mul_pd(x1px2v, EVV)); |
---|
145 | |
---|
146 | _mm_store_pd(&x3[k], ex3v); |
---|
147 | } |
---|
148 | } |
---|
149 | |
---|
150 | scale = 1; |
---|
151 | for(l = 0; scale && (l < 20); l++) |
---|
152 | scale = ((x3[l] < minlikelihood) && (x3[l] > minusminlikelihood)); |
---|
153 | |
---|
154 | if(scale) |
---|
155 | { |
---|
156 | __m128d twoto = _mm_set_pd(twotothe256, twotothe256); |
---|
157 | |
---|
158 | for(l = 0; l < 20; l+=2) |
---|
159 | { |
---|
160 | __m128d ex3v = _mm_mul_pd(_mm_load_pd(&x3[l]),twoto); |
---|
161 | _mm_store_pd(&x3[l], ex3v); |
---|
162 | } |
---|
163 | |
---|
164 | /* |
---|
165 | for(l = 0; l < 20; l++) |
---|
166 | x3[l] *= twotothe256; |
---|
167 | */ |
---|
168 | |
---|
169 | *eVector = *eVector + 1; |
---|
170 | } |
---|
171 | |
---|
172 | return; |
---|
173 | } |
---|
174 | } |
---|
175 | |
---|
176 | static double evaluatePartialGTRCATPROT(int i, double ki, int counter, traversalInfo *ti, double qz, |
---|
177 | int w, double *EIGN, double *EI, double *EV, |
---|
178 | double *tipVector, unsigned char **yVector, |
---|
179 | int branchReference, int mxtips) |
---|
180 | { |
---|
181 | double lz, term; |
---|
182 | double d[20]; |
---|
183 | double *x1, *x2; |
---|
184 | int scale = 0, k, l; |
---|
185 | double |
---|
186 | *lVector = (double *)rax_malloc(sizeof(double) * 20 * mxtips), |
---|
187 | myEI[400] __attribute__ ((aligned (BYTE_ALIGNMENT))); |
---|
188 | |
---|
189 | traversalInfo *trav = &ti[0]; |
---|
190 | |
---|
191 | |
---|
192 | |
---|
193 | for(k = 0; k < 20; k++) |
---|
194 | { |
---|
195 | myEI[k * 20] = 1.0; |
---|
196 | for(l = 1; l < 20; l++) |
---|
197 | myEI[k * 20 + l] = EI[k * 19 + l - 1]; |
---|
198 | } |
---|
199 | |
---|
200 | assert(isTip(trav->pNumber, mxtips)); |
---|
201 | |
---|
202 | x1 = &(tipVector[20 * yVector[trav->pNumber][i]]); |
---|
203 | |
---|
204 | for(k = 1; k < counter; k++) |
---|
205 | computeVectorGTRCATPROT(lVector, &scale, ki, i, ti[k].qz[branchReference], ti[k].rz[branchReference], |
---|
206 | &ti[k], EIGN, myEI, EV, |
---|
207 | tipVector, yVector, mxtips); |
---|
208 | |
---|
209 | x2 = &lVector[20 * (trav->qNumber - mxtips)]; |
---|
210 | |
---|
211 | |
---|
212 | |
---|
213 | assert(0 <= (trav->qNumber - mxtips) && (trav->qNumber - mxtips) < mxtips); |
---|
214 | |
---|
215 | if(qz < zmin) |
---|
216 | lz = zmin; |
---|
217 | lz = log(qz); |
---|
218 | lz *= ki; |
---|
219 | |
---|
220 | d[0] = 1.0; |
---|
221 | for(l = 1; l < 20; l++) |
---|
222 | d[l] = EXP (EIGN[l-1] * lz); |
---|
223 | |
---|
224 | term = 0.0; |
---|
225 | |
---|
226 | for(l = 0; l < 20; l++) |
---|
227 | term += x1[l] * x2[l] * d[l]; |
---|
228 | |
---|
229 | term = LOG(FABS(term)) + (scale * LOG(minlikelihood)); |
---|
230 | |
---|
231 | term = term * w; |
---|
232 | |
---|
233 | rax_free(lVector); |
---|
234 | |
---|
235 | |
---|
236 | return term; |
---|
237 | } |
---|
238 | |
---|
239 | |
---|
240 | #else |
---|
241 | |
---|
242 | static inline void computeVectorGTRCATPROT(double *lVector, int *eVector, double ki, int i, double qz, double rz, |
---|
243 | traversalInfo *ti, double *EIGN, double *EI, double *EV, double *tipVector, |
---|
244 | unsigned char **yVector, int mxtips) |
---|
245 | { |
---|
246 | double *x1, *x2, *x3; |
---|
247 | int |
---|
248 | pNumber = ti->pNumber, |
---|
249 | rNumber = ti->rNumber, |
---|
250 | qNumber = ti->qNumber; |
---|
251 | |
---|
252 | x3 = &(lVector[20 * (pNumber - mxtips)]); |
---|
253 | |
---|
254 | switch(ti->tipCase) |
---|
255 | { |
---|
256 | case TIP_TIP: |
---|
257 | x1 = &(tipVector[20 * yVector[qNumber][i]]); |
---|
258 | x2 = &(tipVector[20 * yVector[rNumber][i]]); |
---|
259 | break; |
---|
260 | case TIP_INNER: |
---|
261 | x1 = &(tipVector[20 * yVector[qNumber][i]]); |
---|
262 | x2 = &( lVector[20 * (rNumber - mxtips)]); |
---|
263 | break; |
---|
264 | case INNER_INNER: |
---|
265 | x1 = &(lVector[20 * (qNumber - mxtips)]); |
---|
266 | x2 = &(lVector[20 * (rNumber - mxtips)]); |
---|
267 | break; |
---|
268 | default: |
---|
269 | assert(0); |
---|
270 | } |
---|
271 | |
---|
272 | { |
---|
273 | double d1[20], d2[20], ump_x1, ump_x2, x1px2, lz1, lz2; |
---|
274 | int l, k, scale; |
---|
275 | |
---|
276 | lz1 = qz * ki; |
---|
277 | lz2 = rz * ki; |
---|
278 | |
---|
279 | for(l = 1; l < 20; l++) |
---|
280 | { |
---|
281 | d1[l] = x1[l] * EXP(EIGN[l - 1] * lz1); |
---|
282 | d2[l] = x2[l] * EXP(EIGN[l - 1] * lz2); |
---|
283 | } |
---|
284 | |
---|
285 | for(l = 0; l < 20; l++) |
---|
286 | x3[l] = 0.0; |
---|
287 | |
---|
288 | for(l = 0; l < 20; l++) |
---|
289 | { |
---|
290 | ump_x1 = x1[0]; |
---|
291 | ump_x2 = x2[0]; |
---|
292 | |
---|
293 | for(k = 1; k < 20; k++) |
---|
294 | { |
---|
295 | ump_x1 += d1[k] * EI[19 * l + k-1]; |
---|
296 | ump_x2 += d2[k] * EI[19 * l + k-1]; |
---|
297 | } |
---|
298 | |
---|
299 | x1px2 = ump_x1 * ump_x2; |
---|
300 | |
---|
301 | for(k = 0; k < 20; k++) |
---|
302 | x3[k] += x1px2 * EV[l * 20 + k]; |
---|
303 | } |
---|
304 | |
---|
305 | scale = 1; |
---|
306 | for(l = 0; scale && (l < 20); l++) |
---|
307 | scale = ((x3[l] < minlikelihood) && (x3[l] > minusminlikelihood)); |
---|
308 | |
---|
309 | if(scale) |
---|
310 | { |
---|
311 | for(l = 0; l < 20; l++) |
---|
312 | x3[l] *= twotothe256; |
---|
313 | *eVector = *eVector + 1; |
---|
314 | } |
---|
315 | |
---|
316 | return; |
---|
317 | } |
---|
318 | } |
---|
319 | |
---|
320 | static double evaluatePartialGTRCATPROT(int i, double ki, int counter, traversalInfo *ti, double qz, |
---|
321 | int w, double *EIGN, double *EI, double *EV, |
---|
322 | double *tipVector, unsigned char **yVector, |
---|
323 | int branchReference, int mxtips) |
---|
324 | { |
---|
325 | double lz, term; |
---|
326 | double d[20]; |
---|
327 | double *x1, *x2; |
---|
328 | int scale = 0, k, l; |
---|
329 | double *lVector = (double *)rax_malloc(sizeof(double) * 20 * mxtips); |
---|
330 | |
---|
331 | traversalInfo *trav = &ti[0]; |
---|
332 | |
---|
333 | assert(isTip(trav->pNumber, mxtips)); |
---|
334 | |
---|
335 | x1 = &(tipVector[20 * yVector[trav->pNumber][i]]); |
---|
336 | |
---|
337 | for(k = 1; k < counter; k++) |
---|
338 | computeVectorGTRCATPROT(lVector, &scale, ki, i, ti[k].qz[branchReference], ti[k].rz[branchReference], |
---|
339 | &ti[k], EIGN, EI, EV, |
---|
340 | tipVector, yVector, mxtips); |
---|
341 | |
---|
342 | x2 = &lVector[20 * (trav->qNumber - mxtips)]; |
---|
343 | |
---|
344 | |
---|
345 | |
---|
346 | assert(0 <= (trav->qNumber - mxtips) && (trav->qNumber - mxtips) < mxtips); |
---|
347 | |
---|
348 | if(qz < zmin) |
---|
349 | lz = zmin; |
---|
350 | lz = log(qz); |
---|
351 | lz *= ki; |
---|
352 | |
---|
353 | d[0] = 1.0; |
---|
354 | for(l = 1; l < 20; l++) |
---|
355 | d[l] = EXP (EIGN[l-1] * lz); |
---|
356 | |
---|
357 | term = 0.0; |
---|
358 | |
---|
359 | for(l = 0; l < 20; l++) |
---|
360 | term += x1[l] * x2[l] * d[l]; |
---|
361 | |
---|
362 | term = LOG(FABS(term)) + (scale * LOG(minlikelihood)); |
---|
363 | |
---|
364 | term = term * w; |
---|
365 | |
---|
366 | rax_free(lVector); |
---|
367 | |
---|
368 | |
---|
369 | return term; |
---|
370 | } |
---|
371 | |
---|
372 | #endif |
---|
373 | |
---|
374 | static inline void computeVectorGTRCATSECONDARY(double *lVector, int *eVector, double ki, int i, double qz, double rz, |
---|
375 | traversalInfo *ti, double *EIGN, double *EI, double *EV, double *tipVector, |
---|
376 | unsigned char **yVector, int mxtips) |
---|
377 | { |
---|
378 | double *x1, *x2, *x3; |
---|
379 | int |
---|
380 | pNumber = ti->pNumber, |
---|
381 | rNumber = ti->rNumber, |
---|
382 | qNumber = ti->qNumber; |
---|
383 | |
---|
384 | x3 = &(lVector[16 * (pNumber - mxtips)]); |
---|
385 | |
---|
386 | switch(ti->tipCase) |
---|
387 | { |
---|
388 | case TIP_TIP: |
---|
389 | x1 = &(tipVector[16 * yVector[qNumber][i]]); |
---|
390 | x2 = &(tipVector[16 * yVector[rNumber][i]]); |
---|
391 | break; |
---|
392 | case TIP_INNER: |
---|
393 | x1 = &(tipVector[16 * yVector[qNumber][i]]); |
---|
394 | x2 = &( lVector[16 * (rNumber - mxtips)]); |
---|
395 | break; |
---|
396 | case INNER_INNER: |
---|
397 | x1 = &(lVector[16 * (qNumber - mxtips)]); |
---|
398 | x2 = &(lVector[16 * (rNumber - mxtips)]); |
---|
399 | break; |
---|
400 | default: |
---|
401 | assert(0); |
---|
402 | } |
---|
403 | |
---|
404 | { |
---|
405 | double d1[16], d2[16], ump_x1, ump_x2, x1px2, lz1, lz2; |
---|
406 | int l, k, scale; |
---|
407 | |
---|
408 | lz1 = qz * ki; |
---|
409 | lz2 = rz * ki; |
---|
410 | |
---|
411 | for(l = 1; l < 16; l++) |
---|
412 | { |
---|
413 | d1[l] = x1[l] * EXP(EIGN[l - 1] * lz1); |
---|
414 | d2[l] = x2[l] * EXP(EIGN[l - 1] * lz2); |
---|
415 | } |
---|
416 | |
---|
417 | for(l = 0; l < 16; l++) |
---|
418 | x3[l] = 0.0; |
---|
419 | |
---|
420 | for(l = 0; l < 16; l++) |
---|
421 | { |
---|
422 | ump_x1 = x1[0]; |
---|
423 | ump_x2 = x2[0]; |
---|
424 | |
---|
425 | for(k = 1; k < 16; k++) |
---|
426 | { |
---|
427 | ump_x1 += d1[k] * EI[15 * l + k-1]; |
---|
428 | ump_x2 += d2[k] * EI[15 * l + k-1]; |
---|
429 | } |
---|
430 | |
---|
431 | x1px2 = ump_x1 * ump_x2; |
---|
432 | |
---|
433 | for(k = 0; k < 16; k++) |
---|
434 | x3[k] += x1px2 * EV[l * 16 + k]; |
---|
435 | } |
---|
436 | |
---|
437 | scale = 1; |
---|
438 | for(l = 0; scale && (l < 16); l++) |
---|
439 | scale = ((x3[l] < minlikelihood) && (x3[l] > minusminlikelihood)); |
---|
440 | |
---|
441 | if(scale) |
---|
442 | { |
---|
443 | for(l = 0; l < 16; l++) |
---|
444 | x3[l] *= twotothe256; |
---|
445 | *eVector = *eVector + 1; |
---|
446 | } |
---|
447 | |
---|
448 | return; |
---|
449 | } |
---|
450 | } |
---|
451 | |
---|
452 | |
---|
453 | static inline void computeVectorFlex(double *lVector, int *eVector, double ki, int i, double qz, double rz, |
---|
454 | traversalInfo *ti, double *EIGN, double *EI, double *EV, double *tipVector, |
---|
455 | unsigned char **yVector, int mxtips, const int numStates) |
---|
456 | { |
---|
457 | double *x1, *x2, *x3; |
---|
458 | int |
---|
459 | pNumber = ti->pNumber, |
---|
460 | rNumber = ti->rNumber, |
---|
461 | qNumber = ti->qNumber; |
---|
462 | |
---|
463 | x3 = &(lVector[numStates * (pNumber - mxtips)]); |
---|
464 | |
---|
465 | switch(ti->tipCase) |
---|
466 | { |
---|
467 | case TIP_TIP: |
---|
468 | x1 = &(tipVector[numStates * yVector[qNumber][i]]); |
---|
469 | x2 = &(tipVector[numStates * yVector[rNumber][i]]); |
---|
470 | break; |
---|
471 | case TIP_INNER: |
---|
472 | x1 = &(tipVector[numStates * yVector[qNumber][i]]); |
---|
473 | x2 = &( lVector[numStates * (rNumber - mxtips)]); |
---|
474 | break; |
---|
475 | case INNER_INNER: |
---|
476 | x1 = &(lVector[numStates * (qNumber - mxtips)]); |
---|
477 | x2 = &(lVector[numStates * (rNumber - mxtips)]); |
---|
478 | break; |
---|
479 | default: |
---|
480 | assert(0); |
---|
481 | } |
---|
482 | |
---|
483 | { |
---|
484 | double d1[64], d2[64], ump_x1, ump_x2, x1px2, lz1, lz2; |
---|
485 | int l, k, scale; |
---|
486 | const int rates = numStates - 1; |
---|
487 | |
---|
488 | lz1 = qz * ki; |
---|
489 | lz2 = rz * ki; |
---|
490 | |
---|
491 | for(l = 1; l < numStates; l++) |
---|
492 | { |
---|
493 | d1[l] = x1[l] * EXP(EIGN[l - 1] * lz1); |
---|
494 | d2[l] = x2[l] * EXP(EIGN[l - 1] * lz2); |
---|
495 | } |
---|
496 | |
---|
497 | for(l = 0; l < numStates; l++) |
---|
498 | x3[l] = 0.0; |
---|
499 | |
---|
500 | for(l = 0; l < numStates; l++) |
---|
501 | { |
---|
502 | ump_x1 = x1[0]; |
---|
503 | ump_x2 = x2[0]; |
---|
504 | |
---|
505 | for(k = 1; k < numStates; k++) |
---|
506 | { |
---|
507 | ump_x1 += d1[k] * EI[rates * l + k-1]; |
---|
508 | ump_x2 += d2[k] * EI[rates * l + k-1]; |
---|
509 | } |
---|
510 | |
---|
511 | x1px2 = ump_x1 * ump_x2; |
---|
512 | |
---|
513 | for(k = 0; k < numStates; k++) |
---|
514 | x3[k] += x1px2 * EV[l * numStates + k]; |
---|
515 | } |
---|
516 | |
---|
517 | scale = 1; |
---|
518 | for(l = 0; scale && (l < numStates); l++) |
---|
519 | scale = ((x3[l] < minlikelihood) && (x3[l] > minusminlikelihood)); |
---|
520 | |
---|
521 | if(scale) |
---|
522 | { |
---|
523 | for(l = 0; l < numStates; l++) |
---|
524 | x3[l] *= twotothe256; |
---|
525 | *eVector = *eVector + 1; |
---|
526 | } |
---|
527 | |
---|
528 | return; |
---|
529 | } |
---|
530 | } |
---|
531 | |
---|
532 | |
---|
533 | static inline void computeVectorGTRCATSECONDARY_6(double *lVector, int *eVector, double ki, int i, double qz, double rz, |
---|
534 | traversalInfo *ti, double *EIGN, double *EI, double *EV, double *tipVector, |
---|
535 | unsigned char **yVector, int mxtips) |
---|
536 | { |
---|
537 | double *x1, *x2, *x3; |
---|
538 | int |
---|
539 | pNumber = ti->pNumber, |
---|
540 | rNumber = ti->rNumber, |
---|
541 | qNumber = ti->qNumber; |
---|
542 | |
---|
543 | x3 = &(lVector[6 * (pNumber - mxtips)]); |
---|
544 | |
---|
545 | switch(ti->tipCase) |
---|
546 | { |
---|
547 | case TIP_TIP: |
---|
548 | x1 = &(tipVector[6 * yVector[qNumber][i]]); |
---|
549 | x2 = &(tipVector[6 * yVector[rNumber][i]]); |
---|
550 | break; |
---|
551 | case TIP_INNER: |
---|
552 | x1 = &(tipVector[6 * yVector[qNumber][i]]); |
---|
553 | x2 = &( lVector[6 * (rNumber - mxtips)]); |
---|
554 | break; |
---|
555 | case INNER_INNER: |
---|
556 | x1 = &(lVector[6 * (qNumber - mxtips)]); |
---|
557 | x2 = &(lVector[6 * (rNumber - mxtips)]); |
---|
558 | break; |
---|
559 | default: |
---|
560 | assert(0); |
---|
561 | } |
---|
562 | |
---|
563 | { |
---|
564 | double d1[6], d2[6], ump_x1, ump_x2, x1px2, lz1, lz2; |
---|
565 | int l, k, scale; |
---|
566 | |
---|
567 | lz1 = qz * ki; |
---|
568 | lz2 = rz * ki; |
---|
569 | |
---|
570 | for(l = 1; l < 6; l++) |
---|
571 | { |
---|
572 | d1[l] = x1[l] * EXP(EIGN[l - 1] * lz1); |
---|
573 | d2[l] = x2[l] * EXP(EIGN[l - 1] * lz2); |
---|
574 | } |
---|
575 | |
---|
576 | for(l = 0; l < 6; l++) |
---|
577 | x3[l] = 0.0; |
---|
578 | |
---|
579 | for(l = 0; l < 6; l++) |
---|
580 | { |
---|
581 | ump_x1 = x1[0]; |
---|
582 | ump_x2 = x2[0]; |
---|
583 | |
---|
584 | for(k = 1; k < 6; k++) |
---|
585 | { |
---|
586 | ump_x1 += d1[k] * EI[5 * l + k-1]; |
---|
587 | ump_x2 += d2[k] * EI[5 * l + k-1]; |
---|
588 | } |
---|
589 | |
---|
590 | x1px2 = ump_x1 * ump_x2; |
---|
591 | |
---|
592 | for(k = 0; k < 6; k++) |
---|
593 | x3[k] += x1px2 * EV[l * 6 + k]; |
---|
594 | } |
---|
595 | |
---|
596 | scale = 1; |
---|
597 | for(l = 0; scale && (l < 6); l++) |
---|
598 | scale = ((x3[l] < minlikelihood) && (x3[l] > minusminlikelihood)); |
---|
599 | |
---|
600 | if(scale) |
---|
601 | { |
---|
602 | for(l = 0; l < 6; l++) |
---|
603 | x3[l] *= twotothe256; |
---|
604 | *eVector = *eVector + 1; |
---|
605 | } |
---|
606 | |
---|
607 | return; |
---|
608 | } |
---|
609 | } |
---|
610 | |
---|
611 | |
---|
612 | static inline void computeVectorGTRCATSECONDARY_7(double *lVector, int *eVector, double ki, int i, double qz, double rz, |
---|
613 | traversalInfo *ti, double *EIGN, double *EI, double *EV, double *tipVector, |
---|
614 | unsigned char **yVector, int mxtips) |
---|
615 | { |
---|
616 | double *x1, *x2, *x3; |
---|
617 | int |
---|
618 | pNumber = ti->pNumber, |
---|
619 | rNumber = ti->rNumber, |
---|
620 | qNumber = ti->qNumber; |
---|
621 | |
---|
622 | x3 = &(lVector[7 * (pNumber - mxtips)]); |
---|
623 | |
---|
624 | switch(ti->tipCase) |
---|
625 | { |
---|
626 | case TIP_TIP: |
---|
627 | x1 = &(tipVector[7 * yVector[qNumber][i]]); |
---|
628 | x2 = &(tipVector[7 * yVector[rNumber][i]]); |
---|
629 | break; |
---|
630 | case TIP_INNER: |
---|
631 | x1 = &(tipVector[7 * yVector[qNumber][i]]); |
---|
632 | x2 = &( lVector[7 * (rNumber - mxtips)]); |
---|
633 | break; |
---|
634 | case INNER_INNER: |
---|
635 | x1 = &(lVector[7 * (qNumber - mxtips)]); |
---|
636 | x2 = &(lVector[7 * (rNumber - mxtips)]); |
---|
637 | break; |
---|
638 | default: |
---|
639 | assert(0); |
---|
640 | } |
---|
641 | |
---|
642 | { |
---|
643 | double d1[7], d2[7], ump_x1, ump_x2, x1px2, lz1, lz2; |
---|
644 | int l, k, scale; |
---|
645 | |
---|
646 | lz1 = qz * ki; |
---|
647 | lz2 = rz * ki; |
---|
648 | |
---|
649 | for(l = 1; l < 7; l++) |
---|
650 | { |
---|
651 | d1[l] = x1[l] * EXP(EIGN[l - 1] * lz1); |
---|
652 | d2[l] = x2[l] * EXP(EIGN[l - 1] * lz2); |
---|
653 | } |
---|
654 | |
---|
655 | for(l = 0; l < 7; l++) |
---|
656 | x3[l] = 0.0; |
---|
657 | |
---|
658 | for(l = 0; l < 7; l++) |
---|
659 | { |
---|
660 | ump_x1 = x1[0]; |
---|
661 | ump_x2 = x2[0]; |
---|
662 | |
---|
663 | for(k = 1; k < 7; k++) |
---|
664 | { |
---|
665 | ump_x1 += d1[k] * EI[6 * l + k-1]; |
---|
666 | ump_x2 += d2[k] * EI[6 * l + k-1]; |
---|
667 | } |
---|
668 | |
---|
669 | x1px2 = ump_x1 * ump_x2; |
---|
670 | |
---|
671 | for(k = 0; k < 7; k++) |
---|
672 | x3[k] += x1px2 * EV[l * 7 + k]; |
---|
673 | } |
---|
674 | |
---|
675 | scale = 1; |
---|
676 | for(l = 0; scale && (l < 7); l++) |
---|
677 | scale = ((x3[l] < minlikelihood) && (x3[l] > minusminlikelihood)); |
---|
678 | |
---|
679 | if(scale) |
---|
680 | { |
---|
681 | for(l = 0; l < 7; l++) |
---|
682 | x3[l] *= twotothe256; |
---|
683 | *eVector = *eVector + 1; |
---|
684 | } |
---|
685 | |
---|
686 | return; |
---|
687 | } |
---|
688 | } |
---|
689 | |
---|
690 | |
---|
691 | static inline void computeVectorGTRCAT(double *lVector, int *eVector, double ki, int i, double qz, double rz, |
---|
692 | traversalInfo *ti, double *EIGN, double *EI, double *EV, double *tipVector, |
---|
693 | unsigned char **yVector, int mxtips) |
---|
694 | { |
---|
695 | double d1[3], d2[3], ump_x1, ump_x2, x1px2[4], lz1, lz2; |
---|
696 | double *x1, *x2, *x3; |
---|
697 | int j, k, |
---|
698 | pNumber = ti->pNumber, |
---|
699 | rNumber = ti->rNumber, |
---|
700 | qNumber = ti->qNumber; |
---|
701 | |
---|
702 | x3 = &lVector[4 * (pNumber - mxtips)]; |
---|
703 | |
---|
704 | |
---|
705 | switch(ti->tipCase) |
---|
706 | { |
---|
707 | case TIP_TIP: |
---|
708 | x1 = &(tipVector[4 * yVector[qNumber][i]]); |
---|
709 | x2 = &(tipVector[4 * yVector[rNumber][i]]); |
---|
710 | break; |
---|
711 | case TIP_INNER: |
---|
712 | x1 = &(tipVector[4 * yVector[qNumber][i]]); |
---|
713 | x2 = &lVector[4 * (rNumber - mxtips)]; |
---|
714 | break; |
---|
715 | case INNER_INNER: |
---|
716 | x1 = &lVector[4 * (qNumber - mxtips)]; |
---|
717 | x2 = &lVector[4 * (rNumber - mxtips)]; |
---|
718 | break; |
---|
719 | default: |
---|
720 | assert(0); |
---|
721 | } |
---|
722 | |
---|
723 | lz1 = qz * ki; |
---|
724 | lz2 = rz * ki; |
---|
725 | |
---|
726 | for(j = 0; j < 3; j++) |
---|
727 | { |
---|
728 | d1[j] = |
---|
729 | x1[j + 1] * |
---|
730 | EXP(EIGN[j] * lz1); |
---|
731 | d2[j] = x2[j + 1] * EXP(EIGN[j] * lz2); |
---|
732 | } |
---|
733 | |
---|
734 | |
---|
735 | for(j = 0; j < 4; j++) |
---|
736 | { |
---|
737 | ump_x1 = x1[0]; |
---|
738 | ump_x2 = x2[0]; |
---|
739 | for(k = 0; k < 3; k++) |
---|
740 | { |
---|
741 | ump_x1 += d1[k] * EI[j * 3 + k]; |
---|
742 | ump_x2 += d2[k] * EI[j * 3 + k]; |
---|
743 | } |
---|
744 | x1px2[j] = ump_x1 * ump_x2; |
---|
745 | } |
---|
746 | |
---|
747 | for(j = 0; j < 4; j++) |
---|
748 | x3[j] = 0.0; |
---|
749 | |
---|
750 | for(j = 0; j < 4; j++) |
---|
751 | for(k = 0; k < 4; k++) |
---|
752 | x3[k] += x1px2[j] * EV[4 * j + k]; |
---|
753 | |
---|
754 | |
---|
755 | if (x3[0] < minlikelihood && x3[0] > minusminlikelihood && |
---|
756 | x3[1] < minlikelihood && x3[1] > minusminlikelihood && |
---|
757 | x3[2] < minlikelihood && x3[2] > minusminlikelihood && |
---|
758 | x3[3] < minlikelihood && x3[3] > minusminlikelihood) |
---|
759 | { |
---|
760 | x3[0] *= twotothe256; |
---|
761 | x3[1] *= twotothe256; |
---|
762 | x3[2] *= twotothe256; |
---|
763 | x3[3] *= twotothe256; |
---|
764 | *eVector = *eVector + 1; |
---|
765 | } |
---|
766 | |
---|
767 | return; |
---|
768 | } |
---|
769 | |
---|
770 | |
---|
771 | |
---|
772 | |
---|
773 | |
---|
774 | |
---|
775 | static inline void computeVectorGTRCAT_BINARY(double *lVector, int *eVector, double ki, int i, double qz, double rz, |
---|
776 | traversalInfo *ti, double *EIGN, double *EI, double *EV, double *tipVector, |
---|
777 | unsigned char **yVector, int mxtips) |
---|
778 | { |
---|
779 | double d1, d2, ump_x1, ump_x2, x1px2[2], lz1, lz2; |
---|
780 | double *x1, *x2, *x3; |
---|
781 | int |
---|
782 | j, k, |
---|
783 | pNumber = ti->pNumber, |
---|
784 | rNumber = ti->rNumber, |
---|
785 | qNumber = ti->qNumber; |
---|
786 | |
---|
787 | x3 = &lVector[2 * (pNumber - mxtips)]; |
---|
788 | |
---|
789 | switch(ti->tipCase) |
---|
790 | { |
---|
791 | case TIP_TIP: |
---|
792 | x1 = &(tipVector[2 * yVector[qNumber][i]]); |
---|
793 | x2 = &(tipVector[2 * yVector[rNumber][i]]); |
---|
794 | break; |
---|
795 | case TIP_INNER: |
---|
796 | x1 = &(tipVector[2 * yVector[qNumber][i]]); |
---|
797 | x2 = &lVector[2 * (rNumber - mxtips)]; |
---|
798 | break; |
---|
799 | case INNER_INNER: |
---|
800 | x1 = &lVector[2 * (qNumber - mxtips)]; |
---|
801 | x2 = &lVector[2 * (rNumber - mxtips)]; |
---|
802 | break; |
---|
803 | default: |
---|
804 | assert(0); |
---|
805 | } |
---|
806 | |
---|
807 | lz1 = qz * ki; |
---|
808 | lz2 = rz * ki; |
---|
809 | |
---|
810 | |
---|
811 | d1 = x1[1] * EXP(EIGN[0] * lz1); |
---|
812 | d2 = x2[1] * EXP(EIGN[0] * lz2); |
---|
813 | |
---|
814 | for(j = 0; j < 2; j++) |
---|
815 | { |
---|
816 | ump_x1 = x1[0]; |
---|
817 | ump_x2 = x2[0]; |
---|
818 | |
---|
819 | ump_x1 += d1 * EI[j]; |
---|
820 | ump_x2 += d2 * EI[j]; |
---|
821 | |
---|
822 | x1px2[j] = ump_x1 * ump_x2; |
---|
823 | } |
---|
824 | |
---|
825 | for(j = 0; j < 2; j++) |
---|
826 | x3[j] = 0.0; |
---|
827 | |
---|
828 | for(j = 0; j < 2; j++) |
---|
829 | for(k = 0; k < 2; k++) |
---|
830 | x3[k] += x1px2[j] * EV[2 * j + k]; |
---|
831 | |
---|
832 | |
---|
833 | if (x3[0] < minlikelihood && x3[0] > minusminlikelihood && |
---|
834 | x3[1] < minlikelihood && x3[1] > minusminlikelihood |
---|
835 | ) |
---|
836 | { |
---|
837 | x3[0] *= twotothe256; |
---|
838 | x3[1] *= twotothe256; |
---|
839 | *eVector = *eVector + 1; |
---|
840 | } |
---|
841 | |
---|
842 | return; |
---|
843 | } |
---|
844 | |
---|
845 | static double evaluatePartialGTRCAT_BINARY(int i, double ki, int counter, traversalInfo *ti, double qz, |
---|
846 | int w, double *EIGN, double *EI, double *EV, |
---|
847 | double *tipVector, unsigned char **yVector, |
---|
848 | int branchReference, int mxtips) |
---|
849 | { |
---|
850 | double lz, term; |
---|
851 | double d; |
---|
852 | double *x1, *x2; |
---|
853 | int scale = 0, k; |
---|
854 | double *lVector = (double *)rax_malloc(sizeof(double) * 2 * mxtips); |
---|
855 | traversalInfo *trav = &ti[0]; |
---|
856 | |
---|
857 | assert(isTip(trav->pNumber, mxtips)); |
---|
858 | |
---|
859 | x1 = &(tipVector[2 * yVector[trav->pNumber][i]]); |
---|
860 | |
---|
861 | for(k = 1; k < counter; k++) |
---|
862 | computeVectorGTRCAT_BINARY(lVector, &scale, ki, i, ti[k].qz[branchReference], ti[k].rz[branchReference], &ti[k], |
---|
863 | EIGN, EI, EV, |
---|
864 | tipVector, yVector, mxtips); |
---|
865 | |
---|
866 | x2 = &lVector[2 * (trav->qNumber - mxtips)]; |
---|
867 | |
---|
868 | |
---|
869 | |
---|
870 | assert(0 <= (trav->qNumber - mxtips) && (trav->qNumber - mxtips) < mxtips); |
---|
871 | |
---|
872 | if(qz < zmin) |
---|
873 | lz = zmin; |
---|
874 | lz = log(qz); |
---|
875 | lz *= ki; |
---|
876 | |
---|
877 | d = EXP (EIGN[0] * lz); |
---|
878 | |
---|
879 | term = x1[0] * x2[0]; |
---|
880 | term += x1[1] * x2[1] * d; |
---|
881 | |
---|
882 | term = LOG(FABS(term)) + (scale * LOG(minlikelihood)); |
---|
883 | |
---|
884 | term = term * w; |
---|
885 | |
---|
886 | rax_free(lVector); |
---|
887 | |
---|
888 | |
---|
889 | return term; |
---|
890 | } |
---|
891 | |
---|
892 | |
---|
893 | static double evaluatePartialGTRCAT(int i, double ki, int counter, traversalInfo *ti, double qz, |
---|
894 | int w, double *EIGN, double *EI, double *EV, |
---|
895 | double *tipVector, unsigned char **yVector, |
---|
896 | int branchReference, int mxtips) |
---|
897 | { |
---|
898 | double lz, term; |
---|
899 | double d[3]; |
---|
900 | double *x1, *x2; |
---|
901 | int scale = 0, k; |
---|
902 | double *lVector = (double *)rax_malloc(sizeof(double) * 4 * mxtips); |
---|
903 | |
---|
904 | traversalInfo *trav = &ti[0]; |
---|
905 | |
---|
906 | assert(isTip(trav->pNumber, mxtips)); |
---|
907 | |
---|
908 | x1 = &(tipVector[4 * yVector[trav->pNumber][i]]); |
---|
909 | |
---|
910 | for(k = 1; k < counter; k++) |
---|
911 | computeVectorGTRCAT(lVector, &scale, ki, i, ti[k].qz[branchReference], ti[k].rz[branchReference], &ti[k], |
---|
912 | EIGN, EI, EV, |
---|
913 | tipVector, yVector, mxtips); |
---|
914 | |
---|
915 | x2 = &lVector[4 * (trav->qNumber - mxtips)]; |
---|
916 | |
---|
917 | assert(0 <= (trav->qNumber - mxtips) && (trav->qNumber - mxtips) < mxtips); |
---|
918 | |
---|
919 | if(qz < zmin) |
---|
920 | lz = zmin; |
---|
921 | lz = log(qz); |
---|
922 | lz *= ki; |
---|
923 | |
---|
924 | d[0] = EXP (EIGN[0] * lz); |
---|
925 | d[1] = EXP (EIGN[1] * lz); |
---|
926 | d[2] = EXP (EIGN[2] * lz); |
---|
927 | |
---|
928 | term = x1[0] * x2[0]; |
---|
929 | term += x1[1] * x2[1] * d[0]; |
---|
930 | term += x1[2] * x2[2] * d[1]; |
---|
931 | term += x1[3] * x2[3] * d[2]; |
---|
932 | |
---|
933 | term = LOG(FABS(term)) + (scale * LOG(minlikelihood)); |
---|
934 | |
---|
935 | term = term * w; |
---|
936 | |
---|
937 | rax_free(lVector); |
---|
938 | |
---|
939 | return term; |
---|
940 | } |
---|
941 | |
---|
942 | |
---|
943 | |
---|
944 | |
---|
945 | static double evaluatePartialGTRCATSECONDARY(int i, double ki, int counter, traversalInfo *ti, double qz, |
---|
946 | int w, double *EIGN, double *EI, double *EV, |
---|
947 | double *tipVector, unsigned char **yVector, |
---|
948 | int branchReference, int mxtips) |
---|
949 | { |
---|
950 | double lz, term; |
---|
951 | double d[16]; |
---|
952 | double *x1, *x2; |
---|
953 | int scale = 0, k, l; |
---|
954 | double *lVector = (double *)rax_malloc(sizeof(double) * 16 * mxtips); |
---|
955 | |
---|
956 | traversalInfo *trav = &ti[0]; |
---|
957 | |
---|
958 | assert(isTip(trav->pNumber, mxtips)); |
---|
959 | |
---|
960 | x1 = &(tipVector[16 * yVector[trav->pNumber][i]]); |
---|
961 | |
---|
962 | for(k = 1; k < counter; k++) |
---|
963 | computeVectorGTRCATSECONDARY(lVector, &scale, ki, i, ti[k].qz[branchReference], ti[k].rz[branchReference], |
---|
964 | &ti[k], EIGN, EI, EV, |
---|
965 | tipVector, yVector, mxtips); |
---|
966 | |
---|
967 | x2 = &lVector[16 * (trav->qNumber - mxtips)]; |
---|
968 | |
---|
969 | |
---|
970 | |
---|
971 | assert(0 <= (trav->qNumber - mxtips) && (trav->qNumber - mxtips) < mxtips); |
---|
972 | |
---|
973 | if(qz < zmin) |
---|
974 | lz = zmin; |
---|
975 | lz = log(qz); |
---|
976 | lz *= ki; |
---|
977 | |
---|
978 | d[0] = 1.0; |
---|
979 | for(l = 1; l < 16; l++) |
---|
980 | d[l] = EXP (EIGN[l-1] * lz); |
---|
981 | |
---|
982 | term = 0.0; |
---|
983 | |
---|
984 | for(l = 0; l < 16; l++) |
---|
985 | term += x1[l] * x2[l] * d[l]; |
---|
986 | |
---|
987 | term = LOG(FABS(term)) + (scale * LOG(minlikelihood)); |
---|
988 | |
---|
989 | term = term * w; |
---|
990 | |
---|
991 | rax_free(lVector); |
---|
992 | |
---|
993 | |
---|
994 | return term; |
---|
995 | } |
---|
996 | |
---|
997 | |
---|
998 | static double evaluatePartialFlex(int i, double ki, int counter, traversalInfo *ti, double qz, |
---|
999 | int w, double *EIGN, double *EI, double *EV, |
---|
1000 | double *tipVector, unsigned char **yVector, |
---|
1001 | int branchReference, int mxtips, const int numStates) |
---|
1002 | { |
---|
1003 | double lz, term; |
---|
1004 | double d[64]; |
---|
1005 | double *x1, *x2; |
---|
1006 | int scale = 0, k, l; |
---|
1007 | double *lVector = (double *)rax_malloc(sizeof(double) * numStates * mxtips); |
---|
1008 | |
---|
1009 | traversalInfo *trav = &ti[0]; |
---|
1010 | |
---|
1011 | assert(isTip(trav->pNumber, mxtips)); |
---|
1012 | |
---|
1013 | x1 = &(tipVector[numStates * yVector[trav->pNumber][i]]); |
---|
1014 | |
---|
1015 | for(k = 1; k < counter; k++) |
---|
1016 | computeVectorFlex(lVector, &scale, ki, i, ti[k].qz[branchReference], ti[k].rz[branchReference], |
---|
1017 | &ti[k], EIGN, EI, EV, |
---|
1018 | tipVector, yVector, mxtips, numStates); |
---|
1019 | |
---|
1020 | x2 = &lVector[numStates * (trav->qNumber - mxtips)]; |
---|
1021 | |
---|
1022 | assert(0 <= (trav->qNumber - mxtips) && (trav->qNumber - mxtips) < mxtips); |
---|
1023 | |
---|
1024 | if(qz < zmin) |
---|
1025 | lz = zmin; |
---|
1026 | lz = log(qz); |
---|
1027 | lz *= ki; |
---|
1028 | |
---|
1029 | d[0] = 1.0; |
---|
1030 | for(l = 1; l < numStates; l++) |
---|
1031 | d[l] = EXP (EIGN[l-1] * lz); |
---|
1032 | |
---|
1033 | term = 0.0; |
---|
1034 | |
---|
1035 | for(l = 0; l < numStates; l++) |
---|
1036 | term += x1[l] * x2[l] * d[l]; |
---|
1037 | |
---|
1038 | term = LOG(FABS(term)) + (scale * LOG(minlikelihood)); |
---|
1039 | |
---|
1040 | term = term * w; |
---|
1041 | |
---|
1042 | rax_free(lVector); |
---|
1043 | |
---|
1044 | |
---|
1045 | return term; |
---|
1046 | } |
---|
1047 | |
---|
1048 | |
---|
1049 | static double evaluatePartialGTRCATSECONDARY_6(int i, double ki, int counter, traversalInfo *ti, double qz, |
---|
1050 | int w, double *EIGN, double *EI, double *EV, |
---|
1051 | double *tipVector, unsigned char **yVector, |
---|
1052 | int branchReference, int mxtips) |
---|
1053 | { |
---|
1054 | double lz, term; |
---|
1055 | double d[6]; |
---|
1056 | double *x1, *x2; |
---|
1057 | int scale = 0, k, l; |
---|
1058 | double *lVector = (double *)rax_malloc(sizeof(double) * 6 * mxtips); |
---|
1059 | |
---|
1060 | traversalInfo *trav = &ti[0]; |
---|
1061 | |
---|
1062 | assert(isTip(trav->pNumber, mxtips)); |
---|
1063 | |
---|
1064 | x1 = &(tipVector[6 * yVector[trav->pNumber][i]]); |
---|
1065 | |
---|
1066 | for(k = 1; k < counter; k++) |
---|
1067 | computeVectorGTRCATSECONDARY_6(lVector, &scale, ki, i, ti[k].qz[branchReference], ti[k].rz[branchReference], |
---|
1068 | &ti[k], EIGN, EI, EV, |
---|
1069 | tipVector, yVector, mxtips); |
---|
1070 | |
---|
1071 | x2 = &lVector[6 * (trav->qNumber - mxtips)]; |
---|
1072 | |
---|
1073 | |
---|
1074 | |
---|
1075 | assert(0 <= (trav->qNumber - mxtips) && (trav->qNumber - mxtips) < mxtips); |
---|
1076 | |
---|
1077 | if(qz < zmin) |
---|
1078 | lz = zmin; |
---|
1079 | lz = log(qz); |
---|
1080 | lz *= ki; |
---|
1081 | |
---|
1082 | d[0] = 1.0; |
---|
1083 | for(l = 1; l < 6; l++) |
---|
1084 | d[l] = EXP (EIGN[l-1] * lz); |
---|
1085 | |
---|
1086 | term = 0.0; |
---|
1087 | |
---|
1088 | for(l = 0; l < 6; l++) |
---|
1089 | term += x1[l] * x2[l] * d[l]; |
---|
1090 | |
---|
1091 | term = LOG(FABS(term)) + (scale * LOG(minlikelihood)); |
---|
1092 | |
---|
1093 | term = term * w; |
---|
1094 | |
---|
1095 | rax_free(lVector); |
---|
1096 | |
---|
1097 | |
---|
1098 | return term; |
---|
1099 | } |
---|
1100 | |
---|
1101 | static double evaluatePartialGTRCATSECONDARY_7(int i, double ki, int counter, traversalInfo *ti, double qz, |
---|
1102 | int w, double *EIGN, double *EI, double *EV, |
---|
1103 | double *tipVector, unsigned char **yVector, |
---|
1104 | int branchReference, int mxtips) |
---|
1105 | { |
---|
1106 | double lz, term; |
---|
1107 | double d[7]; |
---|
1108 | double *x1, *x2; |
---|
1109 | int scale = 0, k, l; |
---|
1110 | double *lVector = (double *)rax_malloc(sizeof(double) * 7 * mxtips); |
---|
1111 | |
---|
1112 | traversalInfo *trav = &ti[0]; |
---|
1113 | |
---|
1114 | assert(isTip(trav->pNumber, mxtips)); |
---|
1115 | |
---|
1116 | x1 = &(tipVector[7 * yVector[trav->pNumber][i]]); |
---|
1117 | |
---|
1118 | for(k = 1; k < counter; k++) |
---|
1119 | computeVectorGTRCATSECONDARY_7(lVector, &scale, ki, i, ti[k].qz[branchReference], ti[k].rz[branchReference], |
---|
1120 | &ti[k], EIGN, EI, EV, |
---|
1121 | tipVector, yVector, mxtips); |
---|
1122 | |
---|
1123 | x2 = &lVector[7 * (trav->qNumber - mxtips)]; |
---|
1124 | |
---|
1125 | |
---|
1126 | |
---|
1127 | assert(0 <= (trav->qNumber - mxtips) && (trav->qNumber - mxtips) < mxtips); |
---|
1128 | |
---|
1129 | if(qz < zmin) |
---|
1130 | lz = zmin; |
---|
1131 | lz = log(qz); |
---|
1132 | lz *= ki; |
---|
1133 | |
---|
1134 | d[0] = 1.0; |
---|
1135 | for(l = 1; l < 7; l++) |
---|
1136 | d[l] = EXP (EIGN[l-1] * lz); |
---|
1137 | |
---|
1138 | term = 0.0; |
---|
1139 | |
---|
1140 | for(l = 0; l < 7; l++) |
---|
1141 | term += x1[l] * x2[l] * d[l]; |
---|
1142 | |
---|
1143 | term = LOG(FABS(term)) + (scale * LOG(minlikelihood)); |
---|
1144 | |
---|
1145 | term = term * w; |
---|
1146 | |
---|
1147 | rax_free(lVector); |
---|
1148 | |
---|
1149 | |
---|
1150 | return term; |
---|
1151 | } |
---|
1152 | |
---|
1153 | |
---|
1154 | |
---|
1155 | |
---|
1156 | /*********************************************************************************************/ |
---|
1157 | |
---|
1158 | |
---|
1159 | |
---|
1160 | void computeFullTraversalInfo(nodeptr p, traversalInfo *ti, int *counter, int maxTips, int numBranches) |
---|
1161 | { |
---|
1162 | if(isTip(p->number, maxTips)) |
---|
1163 | return; |
---|
1164 | |
---|
1165 | { |
---|
1166 | int i; |
---|
1167 | nodeptr q = p->next->back; |
---|
1168 | nodeptr r = p->next->next->back; |
---|
1169 | |
---|
1170 | /* set xnode info at this point */ |
---|
1171 | |
---|
1172 | p->x = 1; |
---|
1173 | p->next->x = 0; |
---|
1174 | p->next->next->x = 0; |
---|
1175 | |
---|
1176 | if(isTip(r->number, maxTips) && isTip(q->number, maxTips)) |
---|
1177 | { |
---|
1178 | ti[*counter].tipCase = TIP_TIP; |
---|
1179 | ti[*counter].pNumber = p->number; |
---|
1180 | ti[*counter].qNumber = q->number; |
---|
1181 | ti[*counter].rNumber = r->number; |
---|
1182 | |
---|
1183 | for(i = 0; i < numBranches; i++) |
---|
1184 | { |
---|
1185 | double z; |
---|
1186 | z = q->z[i]; |
---|
1187 | z = (z > zmin) ? log(z) : log(zmin); |
---|
1188 | ti[*counter].qz[i] = z; |
---|
1189 | |
---|
1190 | z = r->z[i]; |
---|
1191 | z = (z > zmin) ? log(z) : log(zmin); |
---|
1192 | ti[*counter].rz[i] = z; |
---|
1193 | } |
---|
1194 | *counter = *counter + 1; |
---|
1195 | } |
---|
1196 | else |
---|
1197 | { |
---|
1198 | if(isTip(r->number, maxTips) || isTip(q->number, maxTips)) |
---|
1199 | { |
---|
1200 | nodeptr tmp; |
---|
1201 | |
---|
1202 | if(isTip(r->number, maxTips)) |
---|
1203 | { |
---|
1204 | tmp = r; |
---|
1205 | r = q; |
---|
1206 | q = tmp; |
---|
1207 | } |
---|
1208 | |
---|
1209 | computeFullTraversalInfo(r, ti, counter, maxTips, numBranches); |
---|
1210 | |
---|
1211 | ti[*counter].tipCase = TIP_INNER; |
---|
1212 | ti[*counter].pNumber = p->number; |
---|
1213 | ti[*counter].qNumber = q->number; |
---|
1214 | ti[*counter].rNumber = r->number; |
---|
1215 | |
---|
1216 | for(i = 0; i < numBranches; i++) |
---|
1217 | { |
---|
1218 | double z; |
---|
1219 | z = q->z[i]; |
---|
1220 | z = (z > zmin) ? log(z) : log(zmin); |
---|
1221 | ti[*counter].qz[i] = z; |
---|
1222 | |
---|
1223 | z = r->z[i]; |
---|
1224 | z = (z > zmin) ? log(z) : log(zmin); |
---|
1225 | ti[*counter].rz[i] = z; |
---|
1226 | } |
---|
1227 | |
---|
1228 | *counter = *counter + 1; |
---|
1229 | } |
---|
1230 | else |
---|
1231 | { |
---|
1232 | computeFullTraversalInfo(q, ti, counter, maxTips, numBranches); |
---|
1233 | computeFullTraversalInfo(r, ti, counter, maxTips, numBranches); |
---|
1234 | |
---|
1235 | ti[*counter].tipCase = INNER_INNER; |
---|
1236 | ti[*counter].pNumber = p->number; |
---|
1237 | ti[*counter].qNumber = q->number; |
---|
1238 | ti[*counter].rNumber = r->number; |
---|
1239 | for(i = 0; i < numBranches; i++) |
---|
1240 | { |
---|
1241 | double z; |
---|
1242 | z = q->z[i]; |
---|
1243 | z = (z > zmin) ? log(z) : log(zmin); |
---|
1244 | ti[*counter].qz[i] = z; |
---|
1245 | |
---|
1246 | z = r->z[i]; |
---|
1247 | z = (z > zmin) ? log(z) : log(zmin); |
---|
1248 | ti[*counter].rz[i] = z; |
---|
1249 | } |
---|
1250 | |
---|
1251 | *counter = *counter + 1; |
---|
1252 | } |
---|
1253 | } |
---|
1254 | } |
---|
1255 | } |
---|
1256 | |
---|
1257 | void determineFullTraversal(nodeptr p, tree *tr) |
---|
1258 | { |
---|
1259 | nodeptr q = p->back; |
---|
1260 | int k; |
---|
1261 | |
---|
1262 | tr->td[0].ti[0].pNumber = p->number; |
---|
1263 | tr->td[0].ti[0].qNumber = q->number; |
---|
1264 | |
---|
1265 | for(k = 0; k < tr->numBranches; k++) |
---|
1266 | tr->td[0].ti[0].qz[k] = q->z[k]; |
---|
1267 | |
---|
1268 | assert(isTip(p->number, tr->mxtips)); |
---|
1269 | |
---|
1270 | tr->td[0].count = 1; |
---|
1271 | computeFullTraversalInfo(q, &(tr->td[0].ti[0]), &(tr->td[0].count), tr->mxtips, tr->numBranches); |
---|
1272 | computeFullTraversalInfo(p, &(tr->td[0].ti[0]), &(tr->td[0].count), tr->mxtips, tr->numBranches); |
---|
1273 | } |
---|
1274 | |
---|
1275 | |
---|
1276 | |
---|
1277 | |
---|
1278 | |
---|
1279 | double evaluatePartialGeneric (tree *tr, int i, double ki, int _model) |
---|
1280 | { |
---|
1281 | double result; |
---|
1282 | int |
---|
1283 | branchReference, |
---|
1284 | states = tr->partitionData[_model].states; |
---|
1285 | |
---|
1286 | #ifdef _USE_PTHREADS |
---|
1287 | int index = i; |
---|
1288 | #else |
---|
1289 | int index = i - tr->partitionData[_model].lower; |
---|
1290 | #endif |
---|
1291 | |
---|
1292 | if(tr->multiBranch) |
---|
1293 | branchReference = _model; |
---|
1294 | else |
---|
1295 | branchReference = 0; |
---|
1296 | |
---|
1297 | assert(tr->rateHetModel == CAT); |
---|
1298 | |
---|
1299 | |
---|
1300 | |
---|
1301 | switch(tr->partitionData[_model].dataType) |
---|
1302 | { |
---|
1303 | case BINARY_DATA: |
---|
1304 | result = evaluatePartialGTRCAT_BINARY(index, ki, tr->td[0].count, tr->td[0].ti, tr->td[0].ti[0].qz[branchReference], |
---|
1305 | tr->partitionData[_model].wgt[index], |
---|
1306 | tr->partitionData[_model].EIGN, |
---|
1307 | tr->partitionData[_model].EI, |
---|
1308 | tr->partitionData[_model].EV, |
---|
1309 | tr->partitionData[_model].tipVector, |
---|
1310 | tr->partitionData[_model].yVector, branchReference, tr->mxtips); |
---|
1311 | break; |
---|
1312 | case DNA_DATA: |
---|
1313 | result = evaluatePartialGTRCAT(index, ki, tr->td[0].count, tr->td[0].ti, tr->td[0].ti[0].qz[branchReference], |
---|
1314 | tr->partitionData[_model].wgt[index], |
---|
1315 | tr->partitionData[_model].EIGN, |
---|
1316 | tr->partitionData[_model].EI, |
---|
1317 | tr->partitionData[_model].EV, |
---|
1318 | tr->partitionData[_model].tipVector, |
---|
1319 | tr->partitionData[_model].yVector, branchReference, tr->mxtips); |
---|
1320 | break; |
---|
1321 | case AA_DATA: |
---|
1322 | result = evaluatePartialGTRCATPROT(index, ki, tr->td[0].count, tr->td[0].ti, tr->td[0].ti[0].qz[branchReference], |
---|
1323 | tr->partitionData[_model].wgt[index], |
---|
1324 | tr->partitionData[_model].EIGN, |
---|
1325 | tr->partitionData[_model].EI, |
---|
1326 | tr->partitionData[_model].EV, |
---|
1327 | tr->partitionData[_model].tipVector, |
---|
1328 | tr->partitionData[_model].yVector, branchReference, tr->mxtips); |
---|
1329 | break; |
---|
1330 | case SECONDARY_DATA: |
---|
1331 | result = evaluatePartialGTRCATSECONDARY(index, ki, tr->td[0].count, tr->td[0].ti, tr->td[0].ti[0].qz[branchReference], |
---|
1332 | tr->partitionData[_model].wgt[index], |
---|
1333 | tr->partitionData[_model].EIGN, |
---|
1334 | tr->partitionData[_model].EI, |
---|
1335 | tr->partitionData[_model].EV, |
---|
1336 | tr->partitionData[_model].tipVector, |
---|
1337 | tr->partitionData[_model].yVector, branchReference, tr->mxtips); |
---|
1338 | break; |
---|
1339 | case SECONDARY_DATA_6: |
---|
1340 | result = evaluatePartialGTRCATSECONDARY_6(index, ki, tr->td[0].count, tr->td[0].ti, tr->td[0].ti[0].qz[branchReference], |
---|
1341 | tr->partitionData[_model].wgt[index], |
---|
1342 | tr->partitionData[_model].EIGN, |
---|
1343 | tr->partitionData[_model].EI, |
---|
1344 | tr->partitionData[_model].EV, |
---|
1345 | tr->partitionData[_model].tipVector, |
---|
1346 | tr->partitionData[_model].yVector, branchReference, tr->mxtips); |
---|
1347 | break; |
---|
1348 | case SECONDARY_DATA_7: |
---|
1349 | result = evaluatePartialGTRCATSECONDARY_7(index, ki, tr->td[0].count, tr->td[0].ti, tr->td[0].ti[0].qz[branchReference], |
---|
1350 | tr->partitionData[_model].wgt[index], |
---|
1351 | tr->partitionData[_model].EIGN, |
---|
1352 | tr->partitionData[_model].EI, |
---|
1353 | tr->partitionData[_model].EV, |
---|
1354 | tr->partitionData[_model].tipVector, |
---|
1355 | tr->partitionData[_model].yVector, branchReference, tr->mxtips); |
---|
1356 | break; // unintentional fallthrough fixed --ralf |
---|
1357 | case GENERIC_32: |
---|
1358 | result = evaluatePartialFlex(index, ki, tr->td[0].count, tr->td[0].ti, tr->td[0].ti[0].qz[branchReference], |
---|
1359 | tr->partitionData[_model].wgt[index], |
---|
1360 | tr->partitionData[_model].EIGN, |
---|
1361 | tr->partitionData[_model].EI, |
---|
1362 | tr->partitionData[_model].EV, |
---|
1363 | tr->partitionData[_model].tipVector, |
---|
1364 | tr->partitionData[_model].yVector, branchReference, tr->mxtips, states); |
---|
1365 | break; |
---|
1366 | default: |
---|
1367 | assert(0); |
---|
1368 | } |
---|
1369 | |
---|
1370 | |
---|
1371 | return result; |
---|
1372 | } |
---|
1373 | |
---|