Context Navigation

source: branches/port5/GDE/RAxML/evaluatePartialGeneric.c

Visit:

Last change on this file was 5262, checked in by westram, 18 years ago
update to RAxML 7.0.3
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 26.5 KB

Line
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
44
45
46
47	/******************** GTRCAT *************************************/
48
49
50
51
52	#ifdef WIN32
53	static void computeVectorGTRCATPROT(double lVector, int eVector, double ki, int i, double qz, double rz,
54	traversalInfo ti, double EIGN, double EI, double EV, double *tipVector,
55	char **yVector, int mxtips)
56	#else
57	static inline void computeVectorGTRCATPROT(double lVector, int eVector, double ki, int i, double qz, double rz,
58	traversalInfo ti, double EIGN, double EI, double EV, double *tipVector,
59	char **yVector, int mxtips)
60	#endif
61	{
62	double x1, x2, *x3;
63	int ex3,
64	pNumber = ti->pNumber,
65	rNumber = ti->rNumber,
66	qNumber = ti->qNumber;
67
68	x3 = &(lVector[20 * (pNumber - mxtips)]);
69	ex3 = pNumber - mxtips;
70
71	switch(ti->tipCase)
72	{
73	case TIP_TIP:
74	x1 = &(tipVector[20 * yVector[qNumber][i]]);
75	x2 = &(tipVector[20 * yVector[rNumber][i]]);
76	eVector[ex3] = 0;
77	break;
78	case TIP_INNER:
79	x1 = &(tipVector[20 * yVector[qNumber][i]]);
80	x2 = &( lVector[20 * (rNumber - mxtips)]);
81	eVector[ex3] = eVector[rNumber - mxtips];
82	break;
83	case INNER_INNER:
84	x1 = &(lVector[20 * (qNumber - mxtips)]);
85	x2 = &(lVector[20 * (rNumber - mxtips)]);
86	eVector[ex3] = eVector[qNumber - mxtips] + eVector[rNumber - mxtips];
87	break;
88	default:
89	assert(0);
90	}
91
92	{
93	double d1[19], d2[19], ump_x1[20], ump_x2[20], x1px2, lz1, lz2;
94	double left, right, eptr, eptr2;
95	int l, scale;
96
97	lz1 = qz * ki;
98	lz2 = rz * ki;
99
100	left = x1;
101	right = x2;
102
103	d1[0] = left[1] * exp(EIGN[0] * lz1);
104	d1[1] = left[2] * exp(EIGN[1] * lz1);
105	d1[2] = left[3] * exp(EIGN[2] * lz1);
106	d1[3] = left[4] * exp(EIGN[3] * lz1);
107	d1[4] = left[5] * exp(EIGN[4] * lz1);
108	d1[5] = left[6] * exp(EIGN[5] * lz1);
109	d1[6] = left[7] * exp(EIGN[6] * lz1);
110	d1[7] = left[8] * exp(EIGN[7] * lz1);
111	d1[8] = left[9] * exp(EIGN[8] * lz1);
112	d1[9] = left[10] * exp(EIGN[9] * lz1);
113	d1[10] = left[11] * exp(EIGN[10] * lz1);
114	d1[11] = left[12] * exp(EIGN[11] * lz1);
115	d1[12] = left[13] * exp(EIGN[12] * lz1);
116	d1[13] = left[14] * exp(EIGN[13] * lz1);
117	d1[14] = left[15] * exp(EIGN[14] * lz1);
118	d1[15] = left[16] * exp(EIGN[15] * lz1);
119	d1[16] = left[17] * exp(EIGN[16] * lz1);
120	d1[17] = left[18] * exp(EIGN[17] * lz1);
121	d1[18] = left[19] * exp(EIGN[18] * lz1);
122
123	d2[0] = right[1] * exp(EIGN[0] * lz2);
124	d2[1] = right[2] * exp(EIGN[1] * lz2);
125	d2[2] = right[3] * exp(EIGN[2] * lz2);
126	d2[3] = right[4] * exp(EIGN[3] * lz2);
127	d2[4] = right[5] * exp(EIGN[4] * lz2);
128	d2[5] = right[6] * exp(EIGN[5] * lz2);
129	d2[6] = right[7] * exp(EIGN[6] * lz2);
130	d2[7] = right[8] * exp(EIGN[7] * lz2);
131	d2[8] = right[9] * exp(EIGN[8] * lz2);
132	d2[9] = right[10] * exp(EIGN[9] * lz2);
133	d2[10] = right[11] * exp(EIGN[10] * lz2);
134	d2[11] = right[12] * exp(EIGN[11] * lz2);
135	d2[12] = right[13] * exp(EIGN[12] * lz2);
136	d2[13] = right[14] * exp(EIGN[13] * lz2);
137	d2[14] = right[15] * exp(EIGN[14] * lz2);
138	d2[15] = right[16] * exp(EIGN[15] * lz2);
139	d2[16] = right[17] * exp(EIGN[16] * lz2);
140	d2[17] = right[18] * exp(EIGN[17] * lz2);
141	d2[18] = right[19] * exp(EIGN[18] * lz2);
142
143	eptr = EI;
144	for(l = 0; l < 20; l++)
145	{
146	ump_x1[l] = left[0];
147	ump_x1[l] += d1[0] * *eptr++;
148	ump_x1[l] += d1[1] * *eptr++;
149	ump_x1[l] += d1[2] * *eptr++;
150	ump_x1[l] += d1[3] * *eptr++;
151	ump_x1[l] += d1[4] * *eptr++;
152	ump_x1[l] += d1[5] * *eptr++;
153	ump_x1[l] += d1[6] * *eptr++;
154	ump_x1[l] += d1[7] * *eptr++;
155	ump_x1[l] += d1[8] * *eptr++;
156	ump_x1[l] += d1[9] * *eptr++;
157	ump_x1[l] += d1[10] * *eptr++;
158	ump_x1[l] += d1[11] * *eptr++;
159	ump_x1[l] += d1[12] * *eptr++;
160	ump_x1[l] += d1[13] * *eptr++;
161	ump_x1[l] += d1[14] * *eptr++;
162	ump_x1[l] += d1[15] * *eptr++;
163	ump_x1[l] += d1[16] * *eptr++;
164	ump_x1[l] += d1[17] * *eptr++;
165	ump_x1[l] += d1[18] * *eptr++;
166	}
167
168	eptr = EI;
169	for(l = 0; l < 20; l++)
170	{
171	ump_x2[l] = right[0];
172	ump_x2[l] += d2[0] * *eptr++;
173	ump_x2[l] += d2[1] * *eptr++;
174	ump_x2[l] += d2[2] * *eptr++;
175	ump_x2[l] += d2[3] * *eptr++;
176	ump_x2[l] += d2[4] * *eptr++;
177	ump_x2[l] += d2[5] * *eptr++;
178	ump_x2[l] += d2[6] * *eptr++;
179	ump_x2[l] += d2[7] * *eptr++;
180	ump_x2[l] += d2[8] * *eptr++;
181	ump_x2[l] += d2[9] * *eptr++;
182	ump_x2[l] += d2[10] * *eptr++;
183	ump_x2[l] += d2[11] * *eptr++;
184	ump_x2[l] += d2[12] * *eptr++;
185	ump_x2[l] += d2[13] * *eptr++;
186	ump_x2[l] += d2[14] * *eptr++;
187	ump_x2[l] += d2[15] * *eptr++;
188	ump_x2[l] += d2[16] * *eptr++;
189	ump_x2[l] += d2[17] * *eptr++;
190	ump_x2[l] += d2[18] * *eptr++;
191	}
192
193	left = x3;
194	eptr2 = EV;
195
196	x1px2 = ump_x1[0] * ump_x2[0];
197	left[0] = x1px2 * *eptr2++;
198	left[1] = x1px2 * *eptr2++;
199	left[2] = x1px2 * *eptr2++;
200	left[3] = x1px2 * *eptr2++;
201	left[4] = x1px2 * *eptr2++;
202	left[5] = x1px2 * *eptr2++;
203	left[6] = x1px2 * *eptr2++;
204	left[7] = x1px2 * *eptr2++;
205	left[8] = x1px2 * *eptr2++;
206	left[9] = x1px2 * *eptr2++;
207	left[10] = x1px2 * *eptr2++;
208	left[11] = x1px2 * *eptr2++;
209	left[12] = x1px2 * *eptr2++;
210	left[13] = x1px2 * *eptr2++;
211	left[14] = x1px2 * *eptr2++;
212	left[15] = x1px2 * *eptr2++;
213	left[16] = x1px2 * *eptr2++;
214	left[17] = x1px2 * *eptr2++;
215	left[18] = x1px2 * *eptr2++;
216	left[19] = x1px2 * *eptr2++;
217
218	for(l = 1; l < 20; l++)
219	{
220	x1px2 = ump_x1[l] * ump_x2[l];
221	left[0] += x1px2 * *eptr2++;
222	left[1] += x1px2 * *eptr2++;
223	left[2] += x1px2 * *eptr2++;
224	left[3] += x1px2 * *eptr2++;
225	left[4] += x1px2 * *eptr2++;
226	left[5] += x1px2 * *eptr2++;
227	left[6] += x1px2 * *eptr2++;
228	left[7] += x1px2 * *eptr2++;
229	left[8] += x1px2 * *eptr2++;
230	left[9] += x1px2 * *eptr2++;
231	left[10] += x1px2 * *eptr2++;
232	left[11] += x1px2 * *eptr2++;
233	left[12] += x1px2 * *eptr2++;
234	left[13] += x1px2 * *eptr2++;
235	left[14] += x1px2 * *eptr2++;
236	left[15] += x1px2 * *eptr2++;
237	left[16] += x1px2 * *eptr2++;
238	left[17] += x1px2 * *eptr2++;
239	left[18] += x1px2 * *eptr2++;
240	left[19] += x1px2 * *eptr2++;
241	}
242
243	scale = 1;
244	for(l = 0; scale && (l < 20); l++)
245	scale = ((left[l] < minlikelihood) && (left[l] > minusminlikelihood));
246
247	if(scale)
248	{
249	for(l = 0; l < 20; l++)
250	left[l] *= twotothe256;
251	eVector[ex3] += 1;
252	}
253
254	return;
255	}
256	}
257
258	#ifdef WIN32
259	static void computeVectorGTRCAT(double lVector, int eVector, double ki, int i, double qz, double rz,
260	traversalInfo ti, double EIGN, double EI, double EV, double *tipVector,
261	char **yVector, int mxtips)
262	#else
263	static inline void computeVectorGTRCAT(double lVector, int eVector, double ki, int i, double qz, double rz,
264	traversalInfo ti, double EIGN, double EI, double EV, double *tipVector,
265	char **yVector, int mxtips)
266	#endif
267	{
268	double d1[3], d2[3], ump_x1_1, ump_x1_2, ump_x1_3, ump_x1_0,
269	ump_x2_0, ump_x2_1, ump_x2_2, ump_x2_3, x1px2, lz1, lz2;
270	double x1, x2, *x3;
271	int ex3,
272	pNumber = ti->pNumber,
273	rNumber = ti->rNumber,
274	qNumber = ti->qNumber;
275
276	x3 = &lVector[4 * (pNumber - mxtips)];
277	ex3 = pNumber - mxtips;
278
279	switch(ti->tipCase)
280	{
281	case TIP_TIP:
282	x1 = &(tipVector[4 * yVector[qNumber][i]]);
283	x2 = &(tipVector[4 * yVector[rNumber][i]]);
284	eVector[ex3] = 0;
285	break;
286	case TIP_INNER:
287	x1 = &(tipVector[4 * yVector[qNumber][i]]);
288	x2 = &lVector[4 * (rNumber - mxtips)];
289	eVector[ex3] = eVector[rNumber - mxtips];
290	break;
291	case INNER_INNER:
292	x1 = &lVector[4 * (qNumber - mxtips)];
293	x2 = &lVector[4 * (rNumber - mxtips)];
294	eVector[ex3] = eVector[qNumber - mxtips] + eVector[rNumber - mxtips];
295	break;
296	default:
297	assert(0);
298	}
299
300	lz1 = qz * ki;
301	lz2 = rz * ki;
302
303	d1[0] = x1[1] * exp(EIGN[0] * lz1);
304	d2[0] = x2[1] * exp(EIGN[0] * lz2);
305	d1[1] = x1[2] * exp(EIGN[1] * lz1);
306	d2[1] = x2[2] * exp(EIGN[1] * lz2);
307	d1[2] = x1[3] * exp(EIGN[2] * lz1);
308	d2[2] = x2[3] * exp(EIGN[2] * lz2);
309
310	ump_x1_0 = d1[0] * EI[0];
311	ump_x1_0 += d1[1] * EI[1];
312	ump_x1_0 += d1[2] * EI[2];
313	ump_x1_0 += x1[0];
314
315	ump_x1_1 = d1[0] * EI[3];
316	ump_x1_1 += d1[1] * EI[4];
317	ump_x1_1 += d1[2] * EI[5];
318	ump_x1_1 += x1[0];
319
320	ump_x1_2 = d1[0] * EI[6];
321	ump_x1_2 += d1[1] * EI[7];
322	ump_x1_2 += d1[2] * EI[8];
323	ump_x1_2 += x1[0];
324
325	ump_x1_3 = d1[0] * EI[9];
326	ump_x1_3 += d1[1] * EI[10];
327	ump_x1_3 += d1[2] * EI[11];
328	ump_x1_3 += x1[0];
329
330	ump_x2_0 = d2[0] * EI[0];
331	ump_x2_0 += d2[1] * EI[1];
332	ump_x2_0 += d2[2] * EI[2];
333	ump_x2_0 += x2[0];
334
335	ump_x2_1 = d2[0] * EI[3];
336	ump_x2_1 += d2[1] * EI[4];
337	ump_x2_1 += d2[2] * EI[5];
338	ump_x2_1 += x2[0];
339
340	ump_x2_2 = d2[0] * EI[6];
341	ump_x2_2 += d2[1] * EI[7];
342	ump_x2_2 += d2[2] * EI[8];
343	ump_x2_2 += x2[0];
344
345	ump_x2_3 = d2[0] * EI[9];
346	ump_x2_3 += d2[1] * EI[10];
347	ump_x2_3 += d2[2] * EI[11];
348	ump_x2_3 += x2[0];
349
350	x1px2 = ump_x1_0 * ump_x2_0;
351	x3[0] = x1px2 * EV[0];
352	x3[1] = x1px2 * EV[1];
353	x3[2] = x1px2 * EV[2];
354	x3[3] = x1px2 * EV[3];
355
356	x1px2 = ump_x1_1 * ump_x2_1;
357	x3[0] += x1px2 * EV[4];
358	x3[1] += x1px2 * EV[5];
359	x3[2] += x1px2 * EV[6];
360	x3[3] += x1px2 * EV[7];
361
362	x1px2 = ump_x1_2 * ump_x2_2;
363	x3[0] += x1px2 * EV[8];
364	x3[1] += x1px2 * EV[9];
365	x3[2] += x1px2 * EV[10];
366	x3[3] += x1px2 * EV[11];
367
368	x1px2 = ump_x1_3 * ump_x2_3;
369	x3[0] += x1px2 * EV[12];
370	x3[1] += x1px2 * EV[13];
371	x3[2] += x1px2 * EV[14];
372	x3[3] += x1px2 * EV[15];
373
374	if (x3[0] < minlikelihood && x3[0] > minusminlikelihood &&
375	x3[1] < minlikelihood && x3[1] > minusminlikelihood &&
376	x3[2] < minlikelihood && x3[2] > minusminlikelihood &&
377	x3[3] < minlikelihood && x3[3] > minusminlikelihood)
378	{
379	x3[0] *= twotothe256;
380	x3[1] *= twotothe256;
381	x3[2] *= twotothe256;
382	x3[3] *= twotothe256;
383	eVector[ex3] += 1;
384	}
385
386	return;
387	}
388
389
390	static double evaluatePartialGTRCAT(int i, double ki, int counter, traversalInfo *ti, double qz,
391	int w, double EIGN, double EI, double *EV,
392	double tipVector, char *yVector,
393	int branchReference, int mxtips)
394	{
395	double lz, term;
396	double d[3];
397	double x1, x2;
398	int scale, k;
399	double lVector = (double )malloc(sizeof(double) * 4 * mxtips);
400	int eVector = (int )malloc(sizeof(int) * mxtips);
401	traversalInfo *trav = &ti[0];
402
403	assert(isTip(trav->pNumber, mxtips));
404
405	x1 = &(tipVector[4 * yVector[trav->pNumber][i]]);
406
407	for(k = 1; k < counter; k++)
408	computeVectorGTRCAT(lVector, eVector, ki, i, ti[k].qz[branchReference], ti[k].rz[branchReference], &ti[k],
409	EIGN, EI, EV,
410	tipVector, yVector, mxtips);
411
412	x2 = &lVector[4 * (trav->qNumber - mxtips)];
413
414	scale = eVector[trav->qNumber - mxtips];
415
416	assert(0 <= (trav->qNumber - mxtips) && (trav->qNumber - mxtips) < mxtips);
417
418	if(qz < zmin)
419	lz = zmin;
420	lz = log(qz);
421	lz *= ki;
422
423	d[0] = exp (EIGN[0] * lz);
424	d[1] = exp (EIGN[1] * lz);
425	d[2] = exp (EIGN[2] * lz);
426
427	term = x1[0] * x2[0];
428	term += x1[1] * x2[1] * d[0];
429	term += x1[2] * x2[2] * d[1];
430	term += x1[3] * x2[3] * d[2];
431
432	term = log(term) + (scale * log(minlikelihood));
433
434	term = term * w;
435
436	free(lVector);
437	free(eVector);
438
439	return term;
440	}
441
442
443
444	static double evaluatePartialGTRCATPROT(int i, double ki, int counter, traversalInfo *ti, double qz,
445	int w, double EIGN, double EI, double *EV,
446	double tipVector, char *yVector,
447	int branchReference, int mxtips)
448	{
449	double lz, term, left, right;
450	double *ds, d[19];
451	double x1, x2;
452	int scale, k;
453	double lVector = (double )malloc(sizeof(double) * 20 * mxtips);
454	int eVector = (int )malloc(sizeof(int) * mxtips);
455	traversalInfo *trav = &ti[0];
456
457	ds = d;
458
459	assert(isTip(trav->pNumber, mxtips));
460
461	x1 = &(tipVector[20 * yVector[trav->pNumber][i]]);
462
463	for(k = 1; k < counter; k++)
464	computeVectorGTRCATPROT(lVector, eVector, ki, i, ti[k].qz[branchReference], ti[k].rz[branchReference],
465	&ti[k], EIGN, EI, EV,
466	tipVector, yVector, mxtips);
467
468	x2 = &lVector[20 * (trav->qNumber - mxtips)];
469
470	scale = eVector[trav->qNumber - mxtips];
471
472	assert(0 <= (trav->qNumber - mxtips) && (trav->qNumber - mxtips) < mxtips);
473
474	if(qz < zmin)
475	lz = zmin;
476	lz = log(qz);
477	lz *= ki;
478
479	d[0] = exp (EIGN[0] * lz);
480	d[1] = exp (EIGN[1] * lz);
481	d[2] = exp (EIGN[2] * lz);
482	d[3] = exp (EIGN[3] * lz);
483	d[4] = exp (EIGN[4] * lz);
484	d[5] = exp (EIGN[5] * lz);
485	d[6] = exp (EIGN[6] * lz);
486	d[7] = exp (EIGN[7] * lz);
487	d[8] = exp (EIGN[8] * lz);
488	d[9] = exp (EIGN[9] * lz);
489	d[10] = exp (EIGN[10] * lz);
490	d[11] = exp (EIGN[11] * lz);
491	d[12] = exp (EIGN[12] * lz);
492	d[13] = exp (EIGN[13] * lz);
493	d[14] = exp (EIGN[14] * lz);
494	d[15] = exp (EIGN[15] * lz);
495	d[16] = exp (EIGN[16] * lz);
496	d[17] = exp (EIGN[17] * lz);
497	d[18] = exp (EIGN[18] * lz);
498
499
500	left = x1;
501	right = x2;
502
503	term = left++ *right++;
504	term += left++ right++ *ds++;
505	term += left++ right++ *ds++;
506	term += left++ right++ *ds++;
507	term += left++ right++ *ds++;
508
509	term += left++ right++ *ds++;
510	term += left++ right++ *ds++;
511	term += left++ right++ *ds++;
512	term += left++ right++ *ds++;
513	term += left++ right++ *ds++;
514
515	term += left++ right++ *ds++;
516	term += left++ right++ *ds++;
517	term += left++ right++ *ds++;
518	term += left++ right++ *ds++;
519	term += left++ right++ *ds++;
520
521	term += left++ right++ *ds++;
522	term += left++ right++ *ds++;
523	term += left++ right++ *ds++;
524	term += left++ right++ *ds++;
525	term += left++ right++ *ds++;
526
527	term = log(term) + (scale * log(minlikelihood));
528
529	term = term * w;
530
531	free(lVector);
532	free(eVector);
533
534	return term;
535	}
536
537
538
539
540
541
542	/*********************************************************************************************/
543
544
545
546	void computeFullTraversalInfo(nodeptr p, traversalInfo ti, int counter, int maxTips, int numBranches)
547	{
548	if(isTip(p->number, maxTips))
549	return;
550
551	{
552	int i;
553	nodeptr q = p->next->back;
554	nodeptr r = p->next->next->back;
555
556	/* set xnode info at this point */
557
558	p->x = 1;
559	p->next->x = 0;
560	p->next->next->x = 0;
561
562	if(isTip(r->number, maxTips) && isTip(q->number, maxTips))
563	{
564	ti[*counter].tipCase = TIP_TIP;
565	ti[*counter].pNumber = p->number;
566	ti[*counter].qNumber = q->number;
567	ti[*counter].rNumber = r->number;
568
569	for(i = 0; i < numBranches; i++)
570	{
571	double z;
572	z = q->z[i];
573	z = (z > zmin) ? log(z) : log(zmin);
574	ti[*counter].qz[i] = z;
575
576	z = r->z[i];
577	z = (z > zmin) ? log(z) : log(zmin);
578	ti[*counter].rz[i] = z;
579	}
580	counter = counter + 1;
581	}
582	else
583	{
584	if(isTip(r->number, maxTips) \|\| isTip(q->number, maxTips))
585	{
586	nodeptr tmp;
587
588	if(isTip(r->number, maxTips))
589	{
590	tmp = r;
591	r = q;
592	q = tmp;
593	}
594
595	computeFullTraversalInfo(r, ti, counter, maxTips, numBranches);
596
597	ti[*counter].tipCase = TIP_INNER;
598	ti[*counter].pNumber = p->number;
599	ti[*counter].qNumber = q->number;
600	ti[*counter].rNumber = r->number;
601
602	for(i = 0; i < numBranches; i++)
603	{
604	double z;
605	z = q->z[i];
606	z = (z > zmin) ? log(z) : log(zmin);
607	ti[*counter].qz[i] = z;
608
609	z = r->z[i];
610	z = (z > zmin) ? log(z) : log(zmin);
611	ti[*counter].rz[i] = z;
612	}
613
614	counter = counter + 1;
615	}
616	else
617	{
618	computeFullTraversalInfo(q, ti, counter, maxTips, numBranches);
619	computeFullTraversalInfo(r, ti, counter, maxTips, numBranches);
620
621	ti[*counter].tipCase = INNER_INNER;
622	ti[*counter].pNumber = p->number;
623	ti[*counter].qNumber = q->number;
624	ti[*counter].rNumber = r->number;
625	for(i = 0; i < numBranches; i++)
626	{
627	double z;
628	z = q->z[i];
629	z = (z > zmin) ? log(z) : log(zmin);
630	ti[*counter].qz[i] = z;
631
632	z = r->z[i];
633	z = (z > zmin) ? log(z) : log(zmin);
634	ti[*counter].rz[i] = z;
635	}
636
637	counter = counter + 1;
638	}
639	}
640	}
641	}
642
643	void determineFullTraversal(nodeptr p, tree *tr)
644	{
645	nodeptr q = p->back;
646	int k;
647
648	tr->td[0].ti[0].pNumber = p->number;
649	tr->td[0].ti[0].qNumber = q->number;
650
651	for(k = 0; k < tr->numBranches; k++)
652	tr->td[0].ti[0].qz[k] = q->z[k];
653
654	assert(isTip(p->number, tr->rdta->numsp));
655
656	tr->td[0].count = 1;
657	computeFullTraversalInfo(q, &(tr->td[0].ti[0]), &(tr->td[0].count), tr->rdta->numsp, tr->numBranches);
658	computeFullTraversalInfo(p, &(tr->td[0].ti[0]), &(tr->td[0].count), tr->rdta->numsp, tr->numBranches);
659
660
661	}
662
663
664	#ifdef _MULTI_GENE
665
666	void computeFullMultiTraversalInfo(nodeptr p, traversalInfo ti, int counter, int maxTips, int model)
667	{
668	if(isTip(p->number, maxTips))
669	return;
670
671	{
672	/* set xnode info at this point */
673
674	/*p->x = 1;
675	p->next->x = 0;
676	p->next->next->x = 0; */
677
678	if(p->backs[model])
679	{
680	nodeptr q = p->next->backs[model];
681	nodeptr r = p->next->next->backs[model];
682	assert(p == p->next->next->next);
683	p->xs[model] = 1;
684	p->next->xs[model] = 0;
685	p->next->next->xs[model] = 0;
686
687	if(isTip(r->number, maxTips) && isTip(q->number, maxTips))
688	{
689	ti[*counter].tipCase = TIP_TIP;
690	ti[*counter].pNumber = p->number;
691	ti[*counter].qNumber = q->number;
692	ti[*counter].rNumber = r->number;
693
694	{
695	double z;
696	z = q->z[model];
697	z = (z > zmin) ? log(z) : log(zmin);
698	ti[*counter].qz[model] = z;
699
700	z = r->z[model];
701	z = (z > zmin) ? log(z) : log(zmin);
702	ti[*counter].rz[model] = z;
703	}
704
705	counter = counter + 1;
706	}
707	else
708	{
709	if(isTip(r->number, maxTips) \|\| isTip(q->number, maxTips))
710	{
711	nodeptr tmp;
712
713	if(isTip(r->number, maxTips))
714	{
715	tmp = r;
716	r = q;
717	q = tmp;
718	}
719
720	computeFullMultiTraversalInfo(r, ti, counter, maxTips, model);
721
722	ti[*counter].tipCase = TIP_INNER;
723	ti[*counter].pNumber = p->number;
724	ti[*counter].qNumber = q->number;
725	ti[*counter].rNumber = r->number;
726
727
728	{
729	double z;
730	z = q->z[model];
731	z = (z > zmin) ? log(z) : log(zmin);
732	ti[*counter].qz[model] = z;
733
734	z = r->z[model];
735	z = (z > zmin) ? log(z) : log(zmin);
736	ti[*counter].rz[model] = z;
737	}
738
739	counter = counter + 1;
740	}
741	else
742	{
743	computeFullMultiTraversalInfo(q, ti, counter, maxTips, model);
744	computeFullMultiTraversalInfo(r, ti, counter, maxTips, model);
745
746	ti[*counter].tipCase = INNER_INNER;
747	ti[*counter].pNumber = p->number;
748	ti[*counter].qNumber = q->number;
749	ti[*counter].rNumber = r->number;
750
751	{
752	double z;
753	z = q->z[model];
754	z = (z > zmin) ? log(z) : log(zmin);
755	ti[*counter].qz[model] = z;
756
757	z = r->z[model];
758	z = (z > zmin) ? log(z) : log(zmin);
759	ti[*counter].rz[model] = z;
760	}
761
762	counter = counter + 1;
763	}
764	}
765	}
766	else
767	{
768	p->xs[model] = 0;
769	p->next->xs[model] = 0;
770	p->next->next->xs[model] = 0;
771	assert(p == p->next->next->next);
772
773	computeFullMultiTraversalInfo(p->next->back, ti, counter, maxTips, model);
774	computeFullMultiTraversalInfo(p->next->next->back, ti, counter, maxTips, model);
775	}
776	}
777	}
778
779
780	void determineFullMultiTraversal(tree *tr)
781	{
782	int model;
783
784	for(model = 0; model < tr->NumberOfModels; model++)
785	{
786	nodeptr p = tr->startVector[model];
787
788	nodeptr q = p->backs[model];
789
790	tr->td[model].ti[0].pNumber = p->number;
791	tr->td[model].ti[0].qNumber = q->number;
792
793	tr->td[model].ti[0].qz[model] = q->z[model];
794
795	assert(isTip(p->number, tr->rdta->numsp));
796
797	tr->td[model].count = 1;
798
799	computeFullMultiTraversalInfo(q, &(tr->td[model].ti[0]), &(tr->td[model].count), tr->rdta->numsp, model);
800	computeFullMultiTraversalInfo(p, &(tr->td[model].ti[0]), &(tr->td[model].count), tr->rdta->numsp, model);
801
802	/printf("%d %d\n", model, tr->td[model].count);/
803	}
804	}
805	#endif
806
807	#ifdef _LOCAL_DATA
808
809	double evaluatePartialGeneric (tree *localTree, int i, double ki)
810	{
811	double result;
812	int model, branchReference;
813
814	model = localTree->strided_model[i];
815
816	if(localTree->multiBranch)
817	branchReference = model;
818	else
819	branchReference = 0;
820
821	if(localTree->mixedData)
822	{
823	assert(0);
824
825	/*switch(tr->partitionData[model].dataType)
826	{
827	case DNA_DATA:
828	result = evaluatePartialGTRCAT(i, ki, tr->td[0].count, tr->td[0].ti, tr->td[0].ti[0].qz[branchReference],
829	tr->cdta->aliaswgt[i],
830	&(tr->EIGN_DNA[model * 3]), &(tr->EI_DNA[model * 12]), &(tr->EV_DNA[model * 16]),
831	&(tr->tipVectorDNA[model * 64]),
832	tr->yVector, branchReference, tr->mxtips);
833	break;
834	case AA_DATA:
835	result = evaluatePartialGTRCATPROT(i, ki, tr->td[0].count, tr->td[0].ti, tr->td[0].ti[0].qz[branchReference],
836	tr->cdta->aliaswgt[i],
837	&(tr->EIGN_AA[model * 19]), &(tr->EI_AA[model * 380]),
838	&(tr->EV_AA[model * 400]),
839	&(tr->tipVectorAA[model * 460]),
840	tr->yVector, branchReference, tr->mxtips);
841	break;
842	default:
843	assert(0);
844	}
845	*/
846	}
847	else
848	{
849
850	switch(localTree->likelihoodFunction)
851	{
852	case GTRCAT:
853	result = evaluatePartialGTRCAT(i, ki, localTree->td[0].count, localTree->td[0].ti, localTree->td[0].ti[0].qz[0],
854	localTree->strided_aliaswgt[i],
855	&(localTree->EIGN_DNA[0]), &(localTree->EI_DNA[0]), &(localTree->EV_DNA[0]),
856	&(localTree->tipVectorDNA[0]),
857	localTree->strided_yVector, 0, localTree->mxtips);
858	break;
859	case GTRCATMULT:
860	result = evaluatePartialGTRCAT(i, ki, localTree->td[0].count, localTree->td[0].ti, localTree->td[0].ti[0].qz[branchReference],
861	localTree->strided_aliaswgt[i],
862	&(localTree->EIGN_DNA[model * 3]), &(localTree->EI_DNA[model * 12]),
863	&(localTree->EV_DNA[model * 16]),
864	&(localTree->tipVectorDNA[model * 64]),
865	localTree->strided_yVector, branchReference, localTree->mxtips);
866	break;
867	case PROTCAT:
868	result = evaluatePartialGTRCATPROT(i, ki, localTree->td[0].count, localTree->td[0].ti,
869	localTree->td[0].ti[0].qz[0], localTree->strided_aliaswgt[i],
870	&(localTree->EIGN_AA[0]), &(localTree->EI_AA[0]), &(localTree->EV_AA[0]),
871	&(localTree->tipVectorAA[0]),
872	localTree->strided_yVector, 0, localTree->mxtips);
873	break;
874	case PROTCATMULT:
875	result = evaluatePartialGTRCATPROT(i, ki, localTree->td[0].count, localTree->td[0].ti, localTree->td[0].ti[0].qz[branchReference],
876	localTree->strided_aliaswgt[i],
877	&(localTree->EIGN_AA[model * 19]), &(localTree->EI_AA[model * 380]),
878	&(localTree->EV_AA[model * 400]),
879	&(localTree->tipVectorAA[model * 460]),
880	localTree->strided_yVector, branchReference, localTree->mxtips);
881	break;
882	default:
883	assert(0);
884	}
885	}
886
887	return result;
888	}
889
890	#else
891
892	double evaluatePartialGeneric (tree *tr, int i, double ki)
893	{
894	double result;
895	int model, branchReference;
896
897	model = tr->model[i];
898
899	if(tr->multiBranch)
900	branchReference = model;
901	else
902	branchReference = 0;
903
904	if(tr->mixedData)
905	{
906
907
908	switch(tr->partitionData[model].dataType)
909	{
910	case DNA_DATA:
911	result = evaluatePartialGTRCAT(i, ki, tr->td[0].count, tr->td[0].ti, tr->td[0].ti[0].qz[branchReference],
912	tr->cdta->aliaswgt[i],
913	&(tr->EIGN_DNA[model * 3]), &(tr->EI_DNA[model * 12]), &(tr->EV_DNA[model * 16]),
914	&(tr->tipVectorDNA[model * 64]),
915	tr->yVector, branchReference, tr->mxtips);
916	break;
917	case AA_DATA:
918	result = evaluatePartialGTRCATPROT(i, ki, tr->td[0].count, tr->td[0].ti, tr->td[0].ti[0].qz[branchReference],
919	tr->cdta->aliaswgt[i],
920	&(tr->EIGN_AA[model * 19]), &(tr->EI_AA[model * 380]),
921	&(tr->EV_AA[model * 400]),
922	&(tr->tipVectorAA[model * 460]),
923	tr->yVector, branchReference, tr->mxtips);
924	break;
925	default:
926	assert(0);
927	}
928	}
929	else
930	{
931
932	switch(tr->likelihoodFunction)
933	{
934	case GTRCAT:
935	result = evaluatePartialGTRCAT(i, ki, tr->td[0].count, tr->td[0].ti, tr->td[0].ti[0].qz[0], tr->cdta->aliaswgt[i],
936	&(tr->EIGN_DNA[0]), &(tr->EI_DNA[0]), &(tr->EV_DNA[0]),
937	&(tr->tipVectorDNA[0]),
938	tr->yVector, 0, tr->mxtips);
939	break;
940	case GTRCATMULT:
941	result = evaluatePartialGTRCAT(i, ki, tr->td[0].count, tr->td[0].ti, tr->td[0].ti[0].qz[branchReference],
942	tr->cdta->aliaswgt[i],
943	&(tr->EIGN_DNA[model * 3]), &(tr->EI_DNA[model * 12]), &(tr->EV_DNA[model * 16]),
944	&(tr->tipVectorDNA[model * 64]),
945	tr->yVector, branchReference, tr->mxtips);
946	break;
947	case PROTCAT:
948	result = evaluatePartialGTRCATPROT(i, ki, tr->td[0].count, tr->td[0].ti, tr->td[0].ti[0].qz[0], tr->cdta->aliaswgt[i],
949	&(tr->EIGN_AA[0]), &(tr->EI_AA[0]), &(tr->EV_AA[0]),
950	&(tr->tipVectorAA[0]),
951	tr->yVector, 0, tr->mxtips);
952	break;
953	case PROTCATMULT:
954	result = evaluatePartialGTRCATPROT(i, ki, tr->td[0].count, tr->td[0].ti, tr->td[0].ti[0].qz[branchReference],
955	tr->cdta->aliaswgt[i],
956	&(tr->EIGN_AA[model * 19]), &(tr->EI_AA[model * 380]), &(tr->EV_AA[model * 400]),
957	&(tr->tipVectorAA[model * 460]),
958	tr->yVector, branchReference, tr->mxtips);
959	break;
960	default:
961	assert(0);
962	}
963	}
964
965	return result;
966	}
967
968	#endif

Note: See TracBrowser for help on using the repository browser.

Download in other formats: