Context Navigation

← Previous Revision
Next Revision →
Blame
Revision Log

source: branches/profile/GDE/RAxML/evaluateGenericSpecial.c

Visit:

Last change on this file was 10409, checked in by aboeckma, 11 years ago
Updated raxml to current version
Property svn:executable set to ``*
File size: 95.1 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	* This program is free software; you may redistribute it and/or modify its
11	* under the terms of the GNU General Public License as published by the Free
12	* Software Foundation; either version 2 of the License, or (at your option)
13	* any later version.
14	*
15	* This program is distributed in the hope that it will be useful, but
16	* WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
17	* or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
18	* for more details.
19	*
20	*
21	* For any other enquiries send an Email to Alexandros Stamatakis
22	* Alexandros.Stamatakis@epfl.ch
23	*
24	* When publishing work that is based on the results from RAxML-VI-HPC please cite:
25	*
26	* Alexandros Stamatakis:"RAxML-VI-HPC: maximum likelihood-based phylogenetic analyses with thousands of taxa and mixed models".
27	* Bioinformatics 2006; doi: 10.1093/bioinformatics/btl446
28	*/
29
30	#ifndef WIN32
31	#include <unistd.h>
32	#endif
33
34	#include <math.h>
35	#include <time.h>
36	#include <stdlib.h>
37	#include <stdio.h>
38	#include <ctype.h>
39	#include <string.h>
40	#include "axml.h"
41
42
43	#ifdef __SIM_SSE3
44	#include <xmmintrin.h>
45	#include <pmmintrin.h>
46	/#include <tmmintrin.h>/
47	#endif
48
49	#ifdef _USE_PTHREADS
50	extern volatile double *reductionBuffer;
51	extern volatile int NumberOfThreads;
52	#endif
53
54
55
56	extern const unsigned int mask32[32];
57
58
59	static void calcDiagptableFlex(double z, int numberOfCategories, double rptr, double EIGN, double *diagptable, const int numStates)
60	{
61	int
62	i,
63	l;
64
65	double
66	lz,
67	lza[64];
68
69	const int
70	rates = numStates - 1;
71
72	assert(numStates <= 64);
73
74	if (z < zmin)
75	lz = log(zmin);
76	else
77	lz = log(z);
78
79	for(l = 0; l < rates; l++)
80	lza[l] = EIGN[l] * lz;
81
82	for(i = 0; i < numberOfCategories; i++)
83	{
84	diagptable[i * numStates] = 1.0;
85
86	for(l = 1; l < numStates; l++)
87	diagptable[i * numStates + l] = EXP(rptr[i] * lza[l - 1]);
88	}
89	}
90
91
92	static void calcDiagptableFlex_LG4(double z, int numberOfCategories, double rptr, double EIGN[4], double *diagptable, const int numStates)
93	{
94	int
95	i,
96	l;
97
98	double
99	lz;
100
101	assert(numStates <= 64);
102
103	if (z < zmin)
104	lz = log(zmin);
105	else
106	lz = log(z);
107
108	for(i = 0; i < numberOfCategories; i++)
109	{
110	diagptable[i * numStates] = 1.0;
111
112	for(l = 1; l < numStates; l++)
113	diagptable[i * numStates + l] = EXP(rptr[i] * EIGN[i][l - 1] * lz);
114	}
115	}
116
117
118	static double evaluateCatFlex(int ex1, int ex2, int cptr, int wptr,
119	double x1, double x2, double *tipVector,
120	unsigned char tipX1, int n, double diagptable_start, double *vector, boolean writeVector, const boolean fastScaling, const int numStates)
121	{
122	double
123	sum = 0.0,
124	term,
125	*diagptable,
126	*left,
127	*right;
128
129	int
130	i,
131	l;
132
133	if(tipX1)
134	{
135	if(writeVector)
136	for (i = 0; i < n; i++)
137	{
138	left = &(tipVector[numStates * tipX1[i]]);
139	right = &(x2[numStates * i]);
140
141	diagptable = &diagptable_start[numStates * cptr[i]];
142
143	for(l = 0, term = 0.0; l < numStates; l++)
144	term += left[l] * right[l] * diagptable[l];
145
146	if(fastScaling)
147	term = LOG(FABS(term));
148	else
149	term = LOG(FABS(term)) + (ex2[i] * LOG(minlikelihood));
150
151	vector[i] = term;
152
153	sum += wptr[i] * term;
154	}
155	else
156	for (i = 0; i < n; i++)
157	{
158	left = &(tipVector[numStates * tipX1[i]]);
159	right = &(x2[numStates * i]);
160
161	diagptable = &diagptable_start[numStates * cptr[i]];
162
163	for(l = 0, term = 0.0; l < numStates; l++)
164	term += left[l] * right[l] * diagptable[l];
165
166	if(fastScaling)
167	term = LOG(FABS(term));
168	else
169	term = LOG(FABS(term)) + (ex2[i] * LOG(minlikelihood));
170
171	sum += wptr[i] * term;
172	}
173	}
174	else
175	{
176	if(writeVector)
177	for (i = 0; i < n; i++)
178	{
179	left = &x1[numStates * i];
180	right = &x2[numStates * i];
181
182	diagptable = &diagptable_start[numStates * cptr[i]];
183
184	for(l = 0, term = 0.0; l < numStates; l++)
185	term += left[l] * right[l] * diagptable[l];
186
187	if(fastScaling)
188	term = LOG(FABS(term));
189	else
190	term = LOG(FABS(term)) + ((ex1[i] + ex2[i]) * LOG(minlikelihood));
191
192	vector[i] = term;
193
194	sum += wptr[i] * term;
195	}
196	else
197	for (i = 0; i < n; i++)
198	{
199	left = &x1[numStates * i];
200	right = &x2[numStates * i];
201
202	diagptable = &diagptable_start[numStates * cptr[i]];
203
204	for(l = 0, term = 0.0; l < numStates; l++)
205	term += left[l] * right[l] * diagptable[l];
206
207	if(fastScaling)
208	term = LOG(FABS(term));
209	else
210	term = LOG(FABS(term)) + ((ex1[i] + ex2[i]) * LOG(minlikelihood));
211
212	sum += wptr[i] * term;
213	}
214	}
215
216
217
218	return sum;
219	}
220
221
222
223
224	static double evaluateGammaFlex(int ex1, int ex2, int *wptr,
225	double x1, double x2,
226	double *tipVector,
227	unsigned char tipX1, int n, double diagptable, double *vector, boolean writeVector, const boolean fastScaling, const int numStates)
228	{
229	double
230	sum = 0.0,
231	term,
232	*left,
233	*right;
234
235	int
236	i,
237	j,
238	l;
239
240	const int
241	gammaStates = numStates * 4;
242
243	if(tipX1)
244	{
245	if(writeVector)
246	for (i = 0; i < n; i++)
247	{
248	left = &(tipVector[numStates * tipX1[i]]);
249
250	for(j = 0, term = 0.0; j < 4; j++)
251	{
252	right = &(x2[gammaStates * i + numStates * j]);
253
254	for(l = 0; l < numStates; l++)
255	term += left[l] * right[l] * diagptable[j * numStates + l];
256	}
257
258	if(fastScaling)
259	term = LOG(0.25 * FABS(term));
260	else
261	term = LOG(0.25 * FABS(term)) + (ex2[i] * LOG(minlikelihood));
262
263	vector[i] = term;
264
265	sum += wptr[i] * term;
266	}
267	else
268	{
269	for (i = 0; i < n; i++)
270	{
271	left = &(tipVector[numStates * tipX1[i]]);
272
273	for(j = 0, term = 0.0; j < 4; j++)
274	{
275	right = &(x2[gammaStates * i + numStates * j]);
276
277	for(l = 0; l < numStates; l++)
278	term += left[l] * right[l] * diagptable[j * numStates + l];
279	}
280
281	if(fastScaling)
282	term = LOG(0.25 * FABS(term));
283	else
284	term = LOG(0.25 * FABS(term)) + (ex2[i] * LOG(minlikelihood));
285
286	sum += wptr[i] * term;
287	}
288	}
289	}
290	else
291	{
292	if(writeVector)
293	for (i = 0; i < n; i++)
294	{
295
296	for(j = 0, term = 0.0; j < 4; j++)
297	{
298	left = &(x1[gammaStates * i + numStates * j]);
299	right = &(x2[gammaStates * i + numStates * j]);
300
301	for(l = 0; l < numStates; l++)
302	term += left[l] * right[l] * diagptable[j * numStates + l];
303	}
304
305	if(fastScaling)
306	term = LOG(0.25 * FABS(term));
307	else
308	term = LOG(0.25 * FABS(term)) + ((ex1[i] + ex2[i])*LOG(minlikelihood));
309
310	vector[i] = term;
311
312	sum += wptr[i] * term;
313	}
314	else
315	for (i = 0; i < n; i++)
316	{
317
318	for(j = 0, term = 0.0; j < 4; j++)
319	{
320	left = &(x1[gammaStates * i + numStates * j]);
321	right = &(x2[gammaStates * i + numStates * j]);
322
323	for(l = 0; l < numStates; l++)
324	term += left[l] * right[l] * diagptable[j * numStates + l];
325	}
326
327	if(fastScaling)
328	term = LOG(0.25 * FABS(term));
329	else
330	term = LOG(0.25 * FABS(term)) + ((ex1[i] + ex2[i])*LOG(minlikelihood));
331
332	sum += wptr[i] * term;
333	}
334	}
335
336	return sum;
337	}
338
339	static double evaluateGammaFlex_LG4(int ex1, int ex2, int *wptr,
340	double x1, double x2,
341	double *tipVector[4],
342	unsigned char tipX1, int n, double diagptable, double vector, boolean writeVector, const boolean fastScaling, const int numStates, double weights)
343	{
344	double
345	sum = 0.0,
346	term,
347	*left,
348	*right;
349
350	int
351	i,
352	j,
353	l;
354
355	const int
356	gammaStates = numStates * 4;
357
358	if(tipX1)
359	{
360	if(writeVector)
361	for (i = 0; i < n; i++)
362	{
363	for(j = 0, term = 0.0; j < 4; j++)
364	{
365	double
366	t = 0.0;
367
368	left = &(tipVector[j][numStates * tipX1[i]]);
369	right = &(x2[gammaStates * i + numStates * j]);
370
371	for(l = 0; l < numStates; l++)
372	t += left[l] * right[l] * diagptable[j * numStates + l];
373
374	term += weights[j] * t;
375	}
376
377	if(fastScaling)
378	term = LOG(FABS(term));
379	else
380	term = LOG(FABS(term)) + (ex2[i] * LOG(minlikelihood));
381
382	vector[i] = term;
383
384	sum += wptr[i] * term;
385	}
386	else
387	{
388	for (i = 0; i < n; i++)
389	{
390	for(j = 0, term = 0.0; j < 4; j++)
391	{
392	double
393	t = 0.0;
394
395	left = &(tipVector[j][numStates * tipX1[i]]);
396	right = &(x2[gammaStates * i + numStates * j]);
397
398	for(l = 0; l < numStates; l++)
399	t += left[l] * right[l] * diagptable[j * numStates + l];
400
401	term += weights[j] * t;
402	}
403
404	if(fastScaling)
405	term = LOG(FABS(term));
406	else
407	term = LOG(FABS(term)) + (ex2[i] * LOG(minlikelihood));
408
409	sum += wptr[i] * term;
410	}
411	}
412	}
413	else
414	{
415	if(writeVector)
416	for (i = 0; i < n; i++)
417	{
418
419	for(j = 0, term = 0.0; j < 4; j++)
420	{
421	left = &(x1[gammaStates * i + numStates * j]);
422	right = &(x2[gammaStates * i + numStates * j]);
423
424	for(l = 0; l < numStates; l++)
425	term += left[l] * right[l] * diagptable[j * numStates + l];
426	}
427
428	if(fastScaling)
429	term = LOG(0.25 * FABS(term));
430	else
431	term = LOG(0.25 * FABS(term)) + ((ex1[i] + ex2[i])*LOG(minlikelihood));
432
433	vector[i] = term;
434
435	sum += wptr[i] * term;
436	}
437	else
438	for (i = 0; i < n; i++)
439	{
440
441	for(j = 0, term = 0.0; j < 4; j++)
442	{
443	left = &(x1[gammaStates * i + numStates * j]);
444	right = &(x2[gammaStates * i + numStates * j]);
445
446	for(l = 0; l < numStates; l++)
447	term += left[l] * right[l] * diagptable[j * numStates + l];
448	}
449
450	if(fastScaling)
451	term = LOG(0.25 * FABS(term));
452	else
453	term = LOG(0.25 * FABS(term)) + ((ex1[i] + ex2[i])*LOG(minlikelihood));
454
455	sum += wptr[i] * term;
456	}
457	}
458
459	return sum;
460	}
461
462
463
464	static double evaluateGammaInvarFlex (int ex1, int ex2, int wptr, int iptr,
465	double x1, double x2,
466	double tipVector,double tFreqs, double invariants,
467	unsigned char tipX1, int n, double diagptable, double *vector, boolean writeVector, const boolean fastScaling,
468	const int numStates)
469	{
470	double
471	sum = 0.0,
472	term,
473	freqs[64],
474	scaler = 0.25 * (1.0 - invariants),
475	*left,
476	*right;
477
478	int
479	i,
480	j,
481	l;
482
483	const int gammaStates = numStates * 4;
484
485	for(i = 0; i < numStates; i++)
486	freqs[i] = tFreqs[i] * invariants;
487
488	if(tipX1)
489	{
490	if(writeVector)
491	for (i = 0; i < n; i++)
492	{
493	left = &(tipVector[numStates * tipX1[i]]);
494
495	for(j = 0, term = 0.0; j < 4; j++)
496	{
497	right = &(x2[gammaStates * i + numStates * j]);
498
499	for(l = 0; l < numStates; l++)
500	term += left[l] * right[l] * diagptable[j * numStates + l];
501	}
502
503	if(iptr[i] < numStates)
504	if(fastScaling)
505	term = LOG(((scaler * FABS(term)) + freqs[iptr[i]]));
506	else
507	term = LOG(((scaler * FABS(term)) + freqs[iptr[i]])) + ex2[i] * LOG(minlikelihood);
508	else
509	if(fastScaling)
510	term = LOG(scaler * FABS(term));
511	else
512	term = LOG(scaler * FABS(term)) + (ex2[i] * LOG(minlikelihood));
513
514	vector[i] = term;
515
516	sum += wptr[i] * term;
517	}
518	else
519	for (i = 0; i < n; i++)
520	{
521	left = &(tipVector[numStates * tipX1[i]]);
522
523	for(j = 0, term = 0.0; j < 4; j++)
524	{
525	right = &(x2[gammaStates * i + numStates * j]);
526
527	for(l = 0; l < numStates; l++)
528	term += left[l] * right[l] * diagptable[j * numStates + l];
529	}
530
531	if(iptr[i] < numStates)
532	if(fastScaling)
533	term = LOG(((scaler * FABS(term)) + freqs[iptr[i]]));
534	else
535	term = LOG(((scaler * FABS(term)) + freqs[iptr[i]])) + ex2[i] * LOG(minlikelihood);
536	else
537	if(fastScaling)
538	term = LOG(scaler * FABS(term));
539	else
540	term = LOG(scaler * FABS(term)) + (ex2[i] * LOG(minlikelihood));
541
542	sum += wptr[i] * term;
543	}
544	}
545	else
546	{
547	if(writeVector)
548	for (i = 0; i < n; i++)
549	{
550	for(j = 0, term = 0.0; j < 4; j++)
551	{
552	left = &(x1[gammaStates * i + numStates * j]);
553	right = &(x2[gammaStates * i + numStates * j]);
554
555	for(l = 0; l < numStates; l++)
556	term += left[l] * right[l] * diagptable[j * numStates + l];
557	}
558
559	if(iptr[i] < numStates)
560	if(fastScaling)
561	term = LOG(((scaler * FABS(term)) + freqs[iptr[i]]));
562	else
563	term = LOG(((scaler * FABS(term)) + freqs[iptr[i]])) + (ex1[i] + ex2[i]) * LOG(minlikelihood);
564	else
565	if(fastScaling)
566	term = LOG(scaler * FABS(term));
567	else
568	term = LOG(scaler * FABS(term)) + (ex1[i] + ex2[i]) * LOG(minlikelihood);
569
570	vector[i] = term;
571
572	sum += wptr[i] * term;
573	}
574	else
575	for (i = 0; i < n; i++)
576	{
577	for(j = 0, term = 0.0; j < 4; j++)
578	{
579	left = &(x1[gammaStates * i + numStates * j]);
580	right = &(x2[gammaStates * i + numStates * j]);
581
582	for(l = 0; l < numStates; l++)
583	term += left[l] * right[l] * diagptable[j * numStates + l];
584	}
585
586	if(iptr[i] < numStates)
587	if(fastScaling)
588	term = LOG(((scaler * FABS(term)) + freqs[iptr[i]]));
589	else
590	term = LOG(((scaler * FABS(term)) + freqs[iptr[i]])) + (ex1[i] + ex2[i]) * LOG(minlikelihood);
591	else
592	if(fastScaling)
593	term = LOG(scaler * FABS(term));
594	else
595	term = LOG(scaler * FABS(term)) + (ex1[i] + ex2[i]) * LOG(minlikelihood);
596
597	sum += wptr[i] * term;
598	}
599	}
600
601	return sum;
602	}
603
604
605
606	void calcDiagptable(double z, int data, int numberOfCategories, double rptr, double EIGN, double *diagptable)
607	{
608	int i, l;
609	double lz;
610
611	if (z < zmin)
612	lz = log(zmin);
613	else
614	lz = log(z);
615
616	switch(data)
617	{
618	case BINARY_DATA:
619	{
620	double lz1;
621	lz1 = EIGN[0] * lz;
622	for(i = 0; i < numberOfCategories; i++)
623	{
624	diagptable[2 * i] = 1.0;
625	diagptable[2 * i + 1] = EXP(rptr[i] * lz1);
626	}
627	}
628	break;
629	case DNA_DATA:
630	{
631	double lz1, lz2, lz3;
632	lz1 = EIGN[0] * lz;
633	lz2 = EIGN[1] * lz;
634	lz3 = EIGN[2] * lz;
635
636	for(i = 0; i < numberOfCategories; i++)
637	{
638	diagptable[4 * i] = 1.0;
639	diagptable[4 * i + 1] = EXP(rptr[i] * lz1);
640	diagptable[4 * i + 2] = EXP(rptr[i] * lz2);
641	diagptable[4 * i + 3] = EXP(rptr[i] * lz3);
642	}
643	}
644	break;
645	case AA_DATA:
646	{
647	double lza[19];
648
649	for(l = 0; l < 19; l++)
650	lza[l] = EIGN[l] * lz;
651
652	for(i = 0; i < numberOfCategories; i++)
653	{
654	diagptable[i * 20] = 1.0;
655
656	for(l = 1; l < 20; l++)
657	diagptable[i * 20 + l] = EXP(rptr[i] * lza[l - 1]);
658	}
659	}
660	break;
661	case SECONDARY_DATA:
662	{
663	double lza[15];
664
665	for(l = 0; l < 15; l++)
666	lza[l] = EIGN[l] * lz;
667
668	for(i = 0; i < numberOfCategories; i++)
669	{
670	diagptable[i * 16] = 1.0;
671
672	for(l = 1; l < 16; l++)
673	diagptable[i * 16 + l] = EXP(rptr[i] * lza[l - 1]);
674	}
675	}
676	break;
677	case SECONDARY_DATA_6:
678	{
679	double lza[5];
680
681	for(l = 0; l < 5; l++)
682	lza[l] = EIGN[l] * lz;
683
684	for(i = 0; i < numberOfCategories; i++)
685	{
686	diagptable[i * 6] = 1.0;
687
688	for(l = 1; l < 6; l++)
689	diagptable[i * 6 + l] = EXP(rptr[i] * lza[l - 1]);
690	}
691	}
692	break;
693	case SECONDARY_DATA_7:
694	{
695	double lza[6];
696
697	for(l = 0; l < 6; l++)
698	lza[l] = EIGN[l] * lz;
699
700	for(i = 0; i < numberOfCategories; i++)
701	{
702	diagptable[i * 7] = 1.0;
703
704	for(l = 1; l < 7; l++)
705	diagptable[i * 7 + l] = EXP(rptr[i] * lza[l - 1]);
706	}
707	}
708	break;
709	default:
710	assert(0);
711	}
712	}
713
714
715	#ifdef __SIM_SSE3
716
717
718
719
720	static double evaluateGTRCATPROT_SAVE (int ex1, int ex2, int cptr, int wptr,
721	double x1, double x2, double *tipVector,
722	unsigned char tipX1, int n, double diagptable_start, const boolean fastScaling,
723	double x1_gapColumn, double x2_gapColumn, unsigned int x1_gap, unsigned int x2_gap)
724	{
725	double
726	sum = 0.0,
727	term,
728	*diagptable,
729	*left,
730	*right,
731	*left_ptr = x1,
732	*right_ptr = x2;
733
734	int
735	i,
736	l;
737
738	if(tipX1)
739	{
740	for (i = 0; i < n; i++)
741	{
742	left = &(tipVector[20 * tipX1[i]]);
743
744	if(isGap(x2_gap, i))
745	right = x2_gapColumn;
746	else
747	{
748	right = right_ptr;
749	right_ptr += 20;
750	}
751
752	diagptable = &diagptable_start[20 * cptr[i]];
753
754	__m128d tv = _mm_setzero_pd();
755
756	for(l = 0; l < 20; l+=2)
757	{
758	__m128d lv = _mm_load_pd(&left[l]);
759	__m128d rv = _mm_load_pd(&right[l]);
760	__m128d mul = _mm_mul_pd(lv, rv);
761	__m128d dv = _mm_load_pd(&diagptable[l]);
762
763	tv = _mm_add_pd(tv, _mm_mul_pd(mul, dv));
764	}
765
766	tv = _mm_hadd_pd(tv, tv);
767	_mm_storel_pd(&term, tv);
768
769	if(fastScaling)
770	term = LOG(term);
771	else
772	term = LOG(term) + (ex2[i] * LOG(minlikelihood));
773
774	sum += wptr[i] * term;
775	}
776	}
777	else
778	{
779
780	for (i = 0; i < n; i++)
781	{
782	if(isGap(x1_gap, i))
783	left = x1_gapColumn;
784	else
785	{
786	left = left_ptr;
787	left_ptr += 20;
788	}
789
790	if(isGap(x2_gap, i))
791	right = x2_gapColumn;
792	else
793	{
794	right = right_ptr;
795	right_ptr += 20;
796	}
797
798	diagptable = &diagptable_start[20 * cptr[i]];
799
800	__m128d tv = _mm_setzero_pd();
801
802	for(l = 0; l < 20; l+=2)
803	{
804	__m128d lv = _mm_load_pd(&left[l]);
805	__m128d rv = _mm_load_pd(&right[l]);
806	__m128d mul = _mm_mul_pd(lv, rv);
807	__m128d dv = _mm_load_pd(&diagptable[l]);
808
809	tv = _mm_add_pd(tv, _mm_mul_pd(mul, dv));
810	}
811
812	tv = _mm_hadd_pd(tv, tv);
813	_mm_storel_pd(&term, tv);
814
815	if(fastScaling)
816	term = LOG(term);
817	else
818	term = LOG(term) + ((ex1[i] + ex2[i]) * LOG(minlikelihood));
819
820	sum += wptr[i] * term;
821	}
822	}
823
824	return sum;
825	}
826
827
828	static double evaluateGTRCAT_SAVE (int ex1, int ex2, int cptr, int wptr,
829	double x1_start, double x2_start, double *tipVector,
830	unsigned char tipX1, int n, double diagptable_start, const boolean fastScaling,
831	double x1_gapColumn, double x2_gapColumn, unsigned int x1_gap, unsigned int x2_gap)
832	{
833	double sum = 0.0, term;
834	int i;
835
836	double *diagptable,
837	*x1,
838	*x2,
839	*x1_ptr = x1_start,
840	*x2_ptr = x2_start;
841
842	if(tipX1)
843	{
844	for (i = 0; i < n; i++)
845	{
846	double t[2] __attribute__ ((aligned (BYTE_ALIGNMENT)));
847	__m128d x1v1, x1v2, x2v1, x2v2, dv1, dv2;
848
849	x1 = &(tipVector[4 * tipX1[i]]);
850
851	if(isGap(x2_gap, i))
852	x2 = x2_gapColumn;
853	else
854	{
855	x2 = x2_ptr;
856	x2_ptr += 4;
857	}
858
859	diagptable = &diagptable_start[4 * cptr[i]];
860
861	x1v1 = _mm_load_pd(&x1[0]);
862	x1v2 = _mm_load_pd(&x1[2]);
863	x2v1 = _mm_load_pd(&x2[0]);
864	x2v2 = _mm_load_pd(&x2[2]);
865	dv1 = _mm_load_pd(&diagptable[0]);
866	dv2 = _mm_load_pd(&diagptable[2]);
867
868	x1v1 = _mm_mul_pd(x1v1, x2v1);
869	x1v1 = _mm_mul_pd(x1v1, dv1);
870
871	x1v2 = _mm_mul_pd(x1v2, x2v2);
872	x1v2 = _mm_mul_pd(x1v2, dv2);
873
874	x1v1 = _mm_add_pd(x1v1, x1v2);
875
876	_mm_store_pd(t, x1v1);
877
878	if(fastScaling)
879	term = LOG(t[0] + t[1]);
880	else
881	term = LOG(t[0] + t[1]) + (ex2[i] * LOG(minlikelihood));
882
883	sum += wptr[i] * term;
884	}
885	}
886	else
887	{
888	for (i = 0; i < n; i++)
889	{
890	double t[2] __attribute__ ((aligned (BYTE_ALIGNMENT)));
891	__m128d x1v1, x1v2, x2v1, x2v2, dv1, dv2;
892
893	if(isGap(x1_gap, i))
894	x1 = x1_gapColumn;
895	else
896	{
897	x1 = x1_ptr;
898	x1_ptr += 4;
899	}
900
901	if(isGap(x2_gap, i))
902	x2 = x2_gapColumn;
903	else
904	{
905	x2 = x2_ptr;
906	x2_ptr += 4;
907	}
908
909	diagptable = &diagptable_start[4 * cptr[i]];
910
911	x1v1 = _mm_load_pd(&x1[0]);
912	x1v2 = _mm_load_pd(&x1[2]);
913	x2v1 = _mm_load_pd(&x2[0]);
914	x2v2 = _mm_load_pd(&x2[2]);
915	dv1 = _mm_load_pd(&diagptable[0]);
916	dv2 = _mm_load_pd(&diagptable[2]);
917
918	x1v1 = _mm_mul_pd(x1v1, x2v1);
919	x1v1 = _mm_mul_pd(x1v1, dv1);
920
921	x1v2 = _mm_mul_pd(x1v2, x2v2);
922	x1v2 = _mm_mul_pd(x1v2, dv2);
923
924	x1v1 = _mm_add_pd(x1v1, x1v2);
925
926	_mm_store_pd(t, x1v1);
927
928	if(fastScaling)
929	term = LOG(t[0] + t[1]);
930	else
931	term = LOG(t[0] + t[1]) + ((ex1[i] + ex2[i]) * LOG(minlikelihood));
932
933	sum += wptr[i] * term;
934	}
935	}
936
937	return sum;
938	}
939
940	#endif
941
942
943
944
945	static double evaluateGTRCATPROT (int ex1, int ex2, int cptr, int wptr,
946	double x1, double x2, double *tipVector,
947	unsigned char tipX1, int n, double diagptable_start, const boolean fastScaling)
948	{
949	double sum = 0.0, term;
950	double diagptable, left, *right;
951	int i, l;
952
953	if(tipX1)
954	{
955	for (i = 0; i < n; i++)
956	{
957	left = &(tipVector[20 * tipX1[i]]);
958	right = &(x2[20 * i]);
959
960	diagptable = &diagptable_start[20 * cptr[i]];
961	#ifdef __SIM_SSE3
962	__m128d tv = _mm_setzero_pd();
963
964	for(l = 0; l < 20; l+=2)
965	{
966	__m128d lv = _mm_load_pd(&left[l]);
967	__m128d rv = _mm_load_pd(&right[l]);
968	__m128d mul = _mm_mul_pd(lv, rv);
969	__m128d dv = _mm_load_pd(&diagptable[l]);
970
971	tv = _mm_add_pd(tv, _mm_mul_pd(mul, dv));
972	}
973
974	tv = _mm_hadd_pd(tv, tv);
975	_mm_storel_pd(&term, tv);
976	#else
977	for(l = 0, term = 0.0; l < 20; l++)
978	term += left[l] * right[l] * diagptable[l];
979	#endif
980	if(fastScaling)
981	term = LOG(FABS(term));
982	else
983	term = LOG(FABS(term)) + (ex2[i] * LOG(minlikelihood));
984
985	sum += wptr[i] * term;
986	}
987	}
988	else
989	{
990
991	for (i = 0; i < n; i++)
992	{
993	left = &x1[20 * i];
994	right = &x2[20 * i];
995
996	diagptable = &diagptable_start[20 * cptr[i]];
997	#ifdef __SIM_SSE3
998	__m128d tv = _mm_setzero_pd();
999
1000	for(l = 0; l < 20; l+=2)
1001	{
1002	__m128d lv = _mm_load_pd(&left[l]);
1003	__m128d rv = _mm_load_pd(&right[l]);
1004	__m128d mul = _mm_mul_pd(lv, rv);
1005	__m128d dv = _mm_load_pd(&diagptable[l]);
1006
1007	tv = _mm_add_pd(tv, _mm_mul_pd(mul, dv));
1008	}
1009
1010	tv = _mm_hadd_pd(tv, tv);
1011	_mm_storel_pd(&term, tv);
1012	#else
1013	for(l = 0, term = 0.0; l < 20; l++)
1014	term += left[l] * right[l] * diagptable[l];
1015	#endif
1016
1017	if(fastScaling)
1018	term = LOG(FABS(term));
1019	else
1020	term = LOG(FABS(term)) + ((ex1[i] + ex2[i]) * LOG(minlikelihood));
1021
1022	sum += wptr[i] * term;
1023	}
1024	}
1025
1026	return sum;
1027	}
1028
1029
1030	static double evaluateGTRCATSECONDARY (int ex1, int ex2, int cptr, int wptr,
1031	double x1, double x2, double *tipVector,
1032	unsigned char tipX1, int n, double diagptable_start, const boolean fastScaling)
1033	{
1034	double sum = 0.0, term;
1035	double diagptable, left, *right;
1036	int i, l;
1037
1038	if(tipX1)
1039	{
1040	for (i = 0; i < n; i++)
1041	{
1042	left = &(tipVector[16 * tipX1[i]]);
1043	right = &(x2[16 * i]);
1044
1045	diagptable = &diagptable_start[16 * cptr[i]];
1046
1047	for(l = 0, term = 0.0; l < 16; l++)
1048	term += left[l] * right[l] * diagptable[l];
1049
1050	if(fastScaling)
1051	term = LOG(FABS(term));
1052	else
1053	term = LOG(FABS(term)) + (ex2[i] * LOG(minlikelihood));
1054
1055	sum += wptr[i] * term;
1056	}
1057	}
1058	else
1059	{
1060
1061	for (i = 0; i < n; i++)
1062	{
1063	left = &x1[16 * i];
1064	right = &x2[16 * i];
1065
1066	diagptable = &diagptable_start[16 * cptr[i]];
1067
1068	for(l = 0, term = 0.0; l < 16; l++)
1069	term += left[l] * right[l] * diagptable[l];
1070
1071	if(fastScaling)
1072	term = LOG(FABS(term));
1073	else
1074	term = LOG(FABS(term)) + ((ex1[i] + ex2[i]) * LOG(minlikelihood));
1075
1076	sum += wptr[i] * term;
1077	}
1078	}
1079
1080	return sum;
1081	}
1082
1083	static double evaluateGTRCATSECONDARY_6 (int ex1, int ex2, int cptr, int wptr,
1084	double x1, double x2, double *tipVector,
1085	unsigned char tipX1, int n, double diagptable_start, const boolean fastScaling)
1086	{
1087	double sum = 0.0, term;
1088	double diagptable, left, *right;
1089	int i, l;
1090
1091	if(tipX1)
1092	{
1093	for (i = 0; i < n; i++)
1094	{
1095	left = &(tipVector[6 * tipX1[i]]);
1096	right = &(x2[6 * i]);
1097
1098	diagptable = &diagptable_start[6 * cptr[i]];
1099
1100	for(l = 0, term = 0.0; l < 6; l++)
1101	term += left[l] * right[l] * diagptable[l];
1102
1103	if(fastScaling)
1104	term = LOG(FABS(term));
1105	else
1106	term = LOG(FABS(term)) + (ex2[i] * LOG(minlikelihood));
1107
1108	sum += wptr[i] * term;
1109	}
1110	}
1111	else
1112	{
1113
1114	for (i = 0; i < n; i++)
1115	{
1116	left = &x1[6 * i];
1117	right = &x2[6 * i];
1118
1119	diagptable = &diagptable_start[6 * cptr[i]];
1120
1121	for(l = 0, term = 0.0; l < 6; l++)
1122	term += left[l] * right[l] * diagptable[l];
1123
1124	if(fastScaling)
1125	term = LOG(FABS(term));
1126	else
1127	term = LOG(FABS(term)) + ((ex1[i] + ex2[i]) * LOG(minlikelihood));
1128
1129	sum += wptr[i] * term;
1130	}
1131	}
1132
1133	return sum;
1134	}
1135
1136	static double evaluateGTRCATSECONDARY_7(int ex1, int ex2, int cptr, int wptr,
1137	double x1, double x2, double *tipVector,
1138	unsigned char tipX1, int n, double diagptable_start, const boolean fastScaling)
1139	{
1140	double sum = 0.0, term;
1141	double diagptable, left, *right;
1142	int i, l;
1143
1144	if(tipX1)
1145	{
1146	for (i = 0; i < n; i++)
1147	{
1148	left = &(tipVector[7 * tipX1[i]]);
1149	right = &(x2[7 * i]);
1150
1151	diagptable = &diagptable_start[7 * cptr[i]];
1152
1153	for(l = 0, term = 0.0; l < 7; l++)
1154	term += left[l] * right[l] * diagptable[l];
1155
1156	if(fastScaling)
1157	term = LOG(FABS(term));
1158	else
1159	term = LOG(FABS(term)) + (ex2[i] * LOG(minlikelihood));
1160
1161	sum += wptr[i] * term;
1162	}
1163	}
1164	else
1165	{
1166
1167	for (i = 0; i < n; i++)
1168	{
1169	left = &x1[7 * i];
1170	right = &x2[7 * i];
1171
1172	diagptable = &diagptable_start[7 * cptr[i]];
1173
1174	for(l = 0, term = 0.0; l < 7; l++)
1175	term += left[l] * right[l] * diagptable[l];
1176
1177	if(fastScaling)
1178	term = LOG(FABS(term));
1179	else
1180	term = LOG(FABS(term)) + ((ex1[i] + ex2[i]) * LOG(minlikelihood));
1181
1182	sum += wptr[i] * term;
1183	}
1184	}
1185
1186	return sum;
1187	}
1188
1189	static double evaluateGTRCAT_BINARY (int ex1, int ex2, int cptr, int wptr,
1190	double x1_start, double x2_start, double *tipVector,
1191	unsigned char tipX1, int n, double diagptable_start, const boolean fastScaling)
1192	{
1193	double sum = 0.0, term;
1194	int i;
1195	#ifndef __SIM_SSE3
1196	int j;
1197	#endif
1198	double diagptable, x1, *x2;
1199
1200	if(tipX1)
1201	{
1202	for (i = 0; i < n; i++)
1203	{
1204	#ifdef __SIM_SSE3
1205	double t[2] __attribute__ ((aligned (BYTE_ALIGNMENT)));
1206	#endif
1207	x1 = &(tipVector[2 * tipX1[i]]);
1208	x2 = &(x2_start[2 * i]);
1209
1210	diagptable = &(diagptable_start[2 * cptr[i]]);
1211
1212	#ifdef __SIM_SSE3
1213	_mm_store_pd(t, _mm_mul_pd(_mm_load_pd(x1), _mm_mul_pd(_mm_load_pd(x2), _mm_load_pd(diagptable))));
1214
1215	if(fastScaling)
1216	term = LOG(FABS(t[0] + t[1]));
1217	else
1218	term = LOG(FABS(t[0] + t[1])) + (ex2[i] * LOG(minlikelihood));
1219	#else
1220	for(j = 0, term = 0.0; j < 2; j++)
1221	term += x1[j] * x2[j] * diagptable[j];
1222
1223	if(fastScaling)
1224	term = LOG(FABS(term));
1225	else
1226	term = LOG(FABS(term)) + (ex2[i] * LOG(minlikelihood));
1227	#endif
1228
1229	sum += wptr[i] * term;
1230	}
1231	}
1232	else
1233	{
1234	for (i = 0; i < n; i++)
1235	{
1236	#ifdef __SIM_SSE3
1237	double t[2] __attribute__ ((aligned (BYTE_ALIGNMENT)));
1238	#endif
1239	x1 = &x1_start[2 * i];
1240	x2 = &x2_start[2 * i];
1241
1242	diagptable = &diagptable_start[2 * cptr[i]];
1243	#ifdef __SIM_SSE3
1244	_mm_store_pd(t, _mm_mul_pd(_mm_load_pd(x1), _mm_mul_pd(_mm_load_pd(x2), _mm_load_pd(diagptable))));
1245
1246	if(fastScaling)
1247	term = LOG(FABS(t[0] + t[1]));
1248	else
1249	term = LOG(FABS(t[0] + t[1])) + ((ex1[i] + ex2[i]) * LOG(minlikelihood));
1250	#else
1251	for(j = 0, term = 0.0; j < 2; j++)
1252	term += x1[j] * x2[j] * diagptable[j];
1253
1254	if(fastScaling)
1255	term = LOG(FABS(term));
1256	else
1257	term = LOG(FABS(term)) + ((ex1[i] + ex2[i]) * LOG(minlikelihood));
1258	#endif
1259
1260	sum += wptr[i] * term;
1261	}
1262	}
1263
1264	return sum;
1265	}
1266
1267
1268	static double evaluateGTRGAMMA_BINARY(int ex1, int ex2, int *wptr,
1269	double x1_start, double x2_start,
1270	double *tipVector,
1271	unsigned char tipX1, const int n, double diagptable, const boolean fastScaling)
1272	{
1273	double sum = 0.0, term;
1274	int i, j;
1275	#ifndef __SIM_SSE3
1276	int k;
1277	#endif
1278	double x1, x2;
1279
1280	if(tipX1)
1281	{
1282	for (i = 0; i < n; i++)
1283	{
1284	#ifdef __SIM_SSE3
1285	double t[2] __attribute__ ((aligned (BYTE_ALIGNMENT)));
1286	__m128d termv, x1v, x2v, dv;
1287	#endif
1288	x1 = &(tipVector[2 * tipX1[i]]);
1289	x2 = &x2_start[8 * i];
1290	#ifdef __SIM_SSE3
1291	termv = _mm_set1_pd(0.0);
1292
1293	for(j = 0; j < 4; j++)
1294	{
1295	x1v = _mm_load_pd(&x1[0]);
1296	x2v = _mm_load_pd(&x2[j * 2]);
1297	dv = _mm_load_pd(&diagptable[j * 2]);
1298
1299	x1v = _mm_mul_pd(x1v, x2v);
1300	x1v = _mm_mul_pd(x1v, dv);
1301
1302	termv = _mm_add_pd(termv, x1v);
1303	}
1304
1305	_mm_store_pd(t, termv);
1306
1307	if(fastScaling)
1308	term = LOG(0.25 * (FABS(t[0] + t[1])));
1309	else
1310	term = LOG(0.25 * (FABS(t[0] + t[1]))) + (ex2[i] * LOG(minlikelihood));
1311	#else
1312	for(j = 0, term = 0.0; j < 4; j++)
1313	for(k = 0; k < 2; k++)
1314	term += x1[k] * x2[j * 2 + k] * diagptable[j * 2 + k];
1315
1316	if(fastScaling)
1317	term = LOG(0.25 * FABS(term));
1318	else
1319	term = LOG(0.25 * FABS(term)) + ex2[i] * LOG(minlikelihood);
1320	#endif
1321
1322	sum += wptr[i] * term;
1323	}
1324	}
1325	else
1326	{
1327	for (i = 0; i < n; i++)
1328	{
1329	#ifdef __SIM_SSE3
1330	double t[2] __attribute__ ((aligned (BYTE_ALIGNMENT)));
1331	__m128d termv, x1v, x2v, dv;
1332	#endif
1333	x1 = &x1_start[8 * i];
1334	x2 = &x2_start[8 * i];
1335
1336	#ifdef __SIM_SSE3
1337	termv = _mm_set1_pd(0.0);
1338
1339	for(j = 0; j < 4; j++)
1340	{
1341	x1v = _mm_load_pd(&x1[j * 2]);
1342	x2v = _mm_load_pd(&x2[j * 2]);
1343	dv = _mm_load_pd(&diagptable[j * 2]);
1344
1345	x1v = _mm_mul_pd(x1v, x2v);
1346	x1v = _mm_mul_pd(x1v, dv);
1347
1348	termv = _mm_add_pd(termv, x1v);
1349	}
1350
1351	_mm_store_pd(t, termv);
1352
1353
1354	if(fastScaling)
1355	term = LOG(0.25 * (FABS(t[0] + t[1])));
1356	else
1357	term = LOG(0.25 * (FABS(t[0] + t[1]))) + ((ex1[i] +ex2[i]) * LOG(minlikelihood));
1358	#else
1359	for(j = 0, term = 0.0; j < 4; j++)
1360	for(k = 0; k < 2; k++)
1361	term += x1[j * 2 + k] * x2[j * 2 + k] * diagptable[j * 2 + k];
1362
1363	if(fastScaling)
1364	term = LOG(0.25 * FABS(term));
1365	else
1366	term = LOG(0.25 * FABS(term)) + (ex1[i] + ex2[i]) * LOG(minlikelihood);
1367	#endif
1368
1369	sum += wptr[i] * term;
1370	}
1371	}
1372
1373	return sum;
1374	}
1375
1376	static double evaluateGTRGAMMAINVAR_BINARY (int ex1, int ex2, int wptr, int iptr,
1377	double x1_start, double x2_start,
1378	double tipVector, double tFreqs, double invariants,
1379	unsigned char tipX1, int n, double diagptable, const boolean fastScaling)
1380	{
1381	int i, j, k;
1382	double x1, x2;
1383	double
1384	freqs[2],
1385	scaler = 0.25 * (1.0 - invariants),
1386	sum = 0.0,
1387	term;
1388
1389	freqs[0] = tFreqs[0] * invariants;
1390	freqs[1] = tFreqs[1] * invariants;
1391
1392	if(tipX1)
1393	{
1394	for (i = 0; i < n; i++)
1395	{
1396	x1 = &(tipVector[2 * tipX1[i]]);
1397	x2 = &x2_start[8 * i];
1398
1399	for(j = 0, term = 0.0; j < 4; j++)
1400	for(k = 0; k < 2; k++)
1401	term += x1[k] * x2[j * 2 + k] * diagptable[j * 2 + k];
1402
1403	if(iptr[i] < 2)
1404	if(fastScaling)
1405	term = LOG(((scaler * FABS(term)) + freqs[iptr[i]]));
1406	else
1407	term = LOG(((scaler * FABS(term)) + freqs[iptr[i]])) + ex2[i] * LOG(minlikelihood);
1408	else
1409	if(fastScaling)
1410	term = LOG(scaler * FABS(term));
1411	else
1412	term = LOG(scaler * FABS(term)) + (ex2[i] * LOG(minlikelihood));
1413
1414	sum += wptr[i] * term;
1415	}
1416	}
1417	else
1418	{
1419
1420	for (i = 0; i < n; i++)
1421	{
1422	x1 = &x1_start[8 * i];
1423	x2 = &x2_start[8 * i];
1424
1425	for(j = 0, term = 0.0; j < 4; j++)
1426	for(k = 0; k < 2; k++)
1427	term += x1[j * 2 + k] * x2[j * 2 + k] * diagptable[j * 2 + k];
1428
1429	if(iptr[i] < 2)
1430	if(fastScaling)
1431	term = LOG(((scaler * FABS(term)) + freqs[iptr[i]]));
1432	else
1433	term = LOG(((scaler * FABS(term)) + freqs[iptr[i]])) + (ex2[i] + ex1[i]) * LOG(minlikelihood);
1434	else
1435	if(fastScaling)
1436	term = LOG(scaler * FABS(term));
1437	else
1438	term = LOG(scaler * FABS(term)) + ((ex1[i] + ex2[i]) * LOG(minlikelihood));
1439
1440	sum += wptr[i] * term;
1441	}
1442	}
1443
1444	return sum;
1445	}
1446
1447
1448	static double evaluateGTRCAT (int ex1, int ex2, int cptr, int wptr,
1449	double x1_start, double x2_start, double *tipVector,
1450	unsigned char tipX1, int n, double diagptable_start, const boolean fastScaling)
1451	{
1452	double sum = 0.0, term;
1453	int i;
1454	#ifndef __SIM_SSE3
1455	int j;
1456	#endif
1457	double diagptable, x1, *x2;
1458
1459	if(tipX1)
1460	{
1461	for (i = 0; i < n; i++)
1462	{
1463	#ifdef __SIM_SSE3
1464	double t[2] __attribute__ ((aligned (BYTE_ALIGNMENT)));
1465	__m128d x1v1, x1v2, x2v1, x2v2, dv1, dv2;
1466	#endif
1467	x1 = &(tipVector[4 * tipX1[i]]);
1468	x2 = &x2_start[4 * i];
1469
1470	diagptable = &diagptable_start[4 * cptr[i]];
1471
1472	#ifdef __SIM_SSE3
1473	x1v1 = _mm_load_pd(&x1[0]);
1474	x1v2 = _mm_load_pd(&x1[2]);
1475	x2v1 = _mm_load_pd(&x2[0]);
1476	x2v2 = _mm_load_pd(&x2[2]);
1477	dv1 = _mm_load_pd(&diagptable[0]);
1478	dv2 = _mm_load_pd(&diagptable[2]);
1479
1480	x1v1 = _mm_mul_pd(x1v1, x2v1);
1481	x1v1 = _mm_mul_pd(x1v1, dv1);
1482
1483	x1v2 = _mm_mul_pd(x1v2, x2v2);
1484	x1v2 = _mm_mul_pd(x1v2, dv2);
1485
1486	x1v1 = _mm_add_pd(x1v1, x1v2);
1487
1488	_mm_store_pd(t, x1v1);
1489
1490	if(fastScaling)
1491	term = LOG(FABS(t[0] + t[1]));
1492	else
1493	term = LOG(FABS(t[0] + t[1])) + (ex2[i] * LOG(minlikelihood));
1494	#else
1495	for(j = 0, term = 0.0; j < 4; j++)
1496	term += x1[j] * x2[j] * diagptable[j];
1497
1498	/*{
1499	double
1500	term[4],
1501	sum = 0.0;
1502
1503	for(j = 0; j < 4; j++)
1504	{
1505	term[j] = ABS(x1[j] * x2[j] * diagptable[j]);
1506	sum += term[j];
1507	}
1508
1509	printf("RRRRRRR %1.80f %1.80f %1.80f %1.80f\n", term[0]/sum, term[1]/sum, term[2]/sum, term[3]/sum);
1510	}*/
1511
1512	if(fastScaling)
1513	term = LOG(FABS(term));
1514	else
1515	term = LOG(FABS(term)) + (ex2[i] * LOG(minlikelihood));
1516	#endif
1517	sum += wptr[i] * term;
1518	}
1519	}
1520	else
1521	{
1522	for (i = 0; i < n; i++)
1523	{
1524	#ifdef __SIM_SSE3
1525	double t[2] __attribute__ ((aligned (BYTE_ALIGNMENT)));
1526	__m128d x1v1, x1v2, x2v1, x2v2, dv1, dv2;
1527	#endif
1528	x1 = &x1_start[4 * i];
1529	x2 = &x2_start[4 * i];
1530
1531	diagptable = &diagptable_start[4 * cptr[i]];
1532
1533	#ifdef __SIM_SSE3
1534	x1v1 = _mm_load_pd(&x1[0]);
1535	x1v2 = _mm_load_pd(&x1[2]);
1536	x2v1 = _mm_load_pd(&x2[0]);
1537	x2v2 = _mm_load_pd(&x2[2]);
1538	dv1 = _mm_load_pd(&diagptable[0]);
1539	dv2 = _mm_load_pd(&diagptable[2]);
1540
1541	x1v1 = _mm_mul_pd(x1v1, x2v1);
1542	x1v1 = _mm_mul_pd(x1v1, dv1);
1543
1544	x1v2 = _mm_mul_pd(x1v2, x2v2);
1545	x1v2 = _mm_mul_pd(x1v2, dv2);
1546
1547	x1v1 = _mm_add_pd(x1v1, x1v2);
1548
1549	_mm_store_pd(t, x1v1);
1550
1551	if(fastScaling)
1552	term = LOG(FABS(t[0] + t[1]));
1553	else
1554	term = LOG(FABS(t[0] + t[1])) + ((ex1[i] + ex2[i]) * LOG(minlikelihood));
1555	#else
1556
1557	for(j = 0, term = 0.0; j < 4; j++)
1558	term += x1[j] * x2[j] * diagptable[j];
1559
1560	if(fastScaling)
1561	term = LOG(FABS(term));
1562	else
1563	term = LOG(FABS(term)) + ((ex1[i] + ex2[i]) * LOG(minlikelihood));
1564	#endif
1565	sum += wptr[i] * term;
1566	}
1567	}
1568
1569	return sum;
1570	}
1571
1572
1573	#ifdef __SIM_SSE3
1574
1575
1576
1577	static double evaluateGTRGAMMA_GAPPED_SAVE(int ex1, int ex2, int *wptr,
1578	double x1_start, double x2_start,
1579	double *tipVector,
1580	unsigned char tipX1, const int n, double diagptable, const boolean fastScaling,
1581	double x1_gapColumn, double x2_gapColumn, unsigned int x1_gap, unsigned int x2_gap)
1582	{
1583	double sum = 0.0, term;
1584	int i, j;
1585	double
1586	*x1,
1587	*x2,
1588	*x1_ptr = x1_start,
1589	*x2_ptr = x2_start;
1590
1591
1592
1593	if(tipX1)
1594	{
1595
1596
1597	for (i = 0; i < n; i++)
1598	{
1599	double t[2] __attribute__ ((aligned (BYTE_ALIGNMENT)));
1600	__m128d termv, x1v, x2v, dv;
1601
1602	x1 = &(tipVector[4 * tipX1[i]]);
1603	if(x2_gap[i / 32] & mask32[i % 32])
1604	x2 = x2_gapColumn;
1605	else
1606	{
1607	x2 = x2_ptr;
1608	x2_ptr += 16;
1609	}
1610
1611
1612	termv = _mm_set1_pd(0.0);
1613
1614	for(j = 0; j < 4; j++)
1615	{
1616	x1v = _mm_load_pd(&x1[0]);
1617	x2v = _mm_load_pd(&x2[j * 4]);
1618	dv = _mm_load_pd(&diagptable[j * 4]);
1619
1620	x1v = _mm_mul_pd(x1v, x2v);
1621	x1v = _mm_mul_pd(x1v, dv);
1622
1623	termv = _mm_add_pd(termv, x1v);
1624
1625	x1v = _mm_load_pd(&x1[2]);
1626	x2v = _mm_load_pd(&x2[j * 4 + 2]);
1627	dv = _mm_load_pd(&diagptable[j * 4 + 2]);
1628
1629	x1v = _mm_mul_pd(x1v, x2v);
1630	x1v = _mm_mul_pd(x1v, dv);
1631
1632	termv = _mm_add_pd(termv, x1v);
1633	}
1634
1635	_mm_store_pd(t, termv);
1636
1637	if(fastScaling)
1638	term = LOG(0.25 * FABS(t[0] + t[1]));
1639	else
1640	term = LOG(0.25 * FABS(t[0] + t[1])) + (ex2[i] * LOG(minlikelihood));
1641
1642	sum += wptr[i] * term;
1643	}
1644	}
1645	else
1646	{
1647
1648	for (i = 0; i < n; i++)
1649	{
1650
1651	double t[2] __attribute__ ((aligned (BYTE_ALIGNMENT)));
1652	__m128d termv, x1v, x2v, dv;
1653
1654	if(x1_gap[i / 32] & mask32[i % 32])
1655	x1 = x1_gapColumn;
1656	else
1657	{
1658	x1 = x1_ptr;
1659	x1_ptr += 16;
1660	}
1661
1662	if(x2_gap[i / 32] & mask32[i % 32])
1663	x2 = x2_gapColumn;
1664	else
1665	{
1666	x2 = x2_ptr;
1667	x2_ptr += 16;
1668	}
1669
1670	termv = _mm_set1_pd(0.0);
1671
1672	for(j = 0; j < 4; j++)
1673	{
1674	x1v = _mm_load_pd(&x1[j * 4]);
1675	x2v = _mm_load_pd(&x2[j * 4]);
1676	dv = _mm_load_pd(&diagptable[j * 4]);
1677
1678	x1v = _mm_mul_pd(x1v, x2v);
1679	x1v = _mm_mul_pd(x1v, dv);
1680
1681	termv = _mm_add_pd(termv, x1v);
1682
1683	x1v = _mm_load_pd(&x1[j * 4 + 2]);
1684	x2v = _mm_load_pd(&x2[j * 4 + 2]);
1685	dv = _mm_load_pd(&diagptable[j * 4 + 2]);
1686
1687	x1v = _mm_mul_pd(x1v, x2v);
1688	x1v = _mm_mul_pd(x1v, dv);
1689
1690	termv = _mm_add_pd(termv, x1v);
1691	}
1692
1693	_mm_store_pd(t, termv);
1694
1695	if(fastScaling)
1696	term = LOG(0.25 * FABS(t[0] + t[1]));
1697	else
1698	term = LOG(0.25 * FABS(t[0] + t[1])) + ((ex1[i] + ex2[i]) * LOG(minlikelihood));
1699
1700	sum += wptr[i] * term;
1701	}
1702	}
1703
1704	return sum;
1705	}
1706
1707	#else
1708
1709
1710
1711	#endif
1712
1713	static double evaluateGTRGAMMA(int ex1, int ex2, int *wptr,
1714	double x1_start, double x2_start,
1715	double *tipVector,
1716	unsigned char tipX1, const int n, double diagptable, const boolean fastScaling)
1717	{
1718	double sum = 0.0, term;
1719	int i, j;
1720	#ifndef __SIM_SSE3
1721	int k;
1722	#endif
1723	double x1, x2;
1724
1725
1726
1727	if(tipX1)
1728	{
1729	for (i = 0; i < n; i++)
1730	{
1731	#ifdef __SIM_SSE3
1732	double t[2] __attribute__ ((aligned (BYTE_ALIGNMENT)));
1733	__m128d termv, x1v, x2v, dv;
1734	#endif
1735	x1 = &(tipVector[4 * tipX1[i]]);
1736	x2 = &x2_start[16 * i];
1737
1738	#ifdef __SIM_SSE3
1739	termv = _mm_set1_pd(0.0);
1740
1741	for(j = 0; j < 4; j++)
1742	{
1743	x1v = _mm_load_pd(&x1[0]);
1744	x2v = _mm_load_pd(&x2[j * 4]);
1745	dv = _mm_load_pd(&diagptable[j * 4]);
1746
1747	x1v = _mm_mul_pd(x1v, x2v);
1748	x1v = _mm_mul_pd(x1v, dv);
1749
1750	termv = _mm_add_pd(termv, x1v);
1751
1752	x1v = _mm_load_pd(&x1[2]);
1753	x2v = _mm_load_pd(&x2[j * 4 + 2]);
1754	dv = _mm_load_pd(&diagptable[j * 4 + 2]);
1755
1756	x1v = _mm_mul_pd(x1v, x2v);
1757	x1v = _mm_mul_pd(x1v, dv);
1758
1759	termv = _mm_add_pd(termv, x1v);
1760	}
1761
1762	_mm_store_pd(t, termv);
1763
1764
1765	if(fastScaling)
1766	term = LOG(0.25 * FABS(t[0] + t[1]));
1767	else
1768	term = LOG(0.25 * FABS(t[0] + t[1])) + (ex2[i] * LOG(minlikelihood));
1769	#else
1770	for(j = 0, term = 0.0; j < 4; j++)
1771	for(k = 0; k < 4; k++)
1772	term += x1[k] * x2[j * 4 + k] * diagptable[j * 4 + k];
1773
1774	if(fastScaling)
1775	term = LOG(0.25 * FABS(term));
1776	else
1777	term = LOG(0.25 * FABS(term)) + ex2[i] * LOG(minlikelihood);
1778	#endif
1779
1780	sum += wptr[i] * term;
1781	}
1782	}
1783	else
1784	{
1785	for (i = 0; i < n; i++)
1786	{
1787	#ifdef __SIM_SSE3
1788	double t[2] __attribute__ ((aligned (BYTE_ALIGNMENT)));
1789	__m128d termv, x1v, x2v, dv;
1790	#endif
1791
1792	x1 = &x1_start[16 * i];
1793	x2 = &x2_start[16 * i];
1794
1795	#ifdef __SIM_SSE3
1796	termv = _mm_set1_pd(0.0);
1797
1798	for(j = 0; j < 4; j++)
1799	{
1800	x1v = _mm_load_pd(&x1[j * 4]);
1801	x2v = _mm_load_pd(&x2[j * 4]);
1802	dv = _mm_load_pd(&diagptable[j * 4]);
1803
1804	x1v = _mm_mul_pd(x1v, x2v);
1805	x1v = _mm_mul_pd(x1v, dv);
1806
1807	termv = _mm_add_pd(termv, x1v);
1808
1809	x1v = _mm_load_pd(&x1[j * 4 + 2]);
1810	x2v = _mm_load_pd(&x2[j * 4 + 2]);
1811	dv = _mm_load_pd(&diagptable[j * 4 + 2]);
1812
1813	x1v = _mm_mul_pd(x1v, x2v);
1814	x1v = _mm_mul_pd(x1v, dv);
1815
1816	termv = _mm_add_pd(termv, x1v);
1817	}
1818
1819	_mm_store_pd(t, termv);
1820
1821	if(fastScaling)
1822	term = LOG(0.25 * FABS(t[0] + t[1]));
1823	else
1824	term = LOG(0.25 * FABS(t[0] + t[1])) + ((ex1[i] + ex2[i]) * LOG(minlikelihood));
1825	#else
1826	for(j = 0, term = 0.0; j < 4; j++)
1827	for(k = 0; k < 4; k++)
1828	term += x1[j * 4 + k] * x2[j * 4 + k] * diagptable[j * 4 + k];
1829
1830	if(fastScaling)
1831	term = LOG(0.25 * FABS(term));
1832	else
1833	term = LOG(0.25 * FABS(term)) + (ex1[i] + ex2[i]) * LOG(minlikelihood);
1834	#endif
1835
1836	sum += wptr[i] * term;
1837	}
1838	}
1839
1840	return sum;
1841	}
1842
1843
1844
1845
1846
1847
1848
1849
1850
1851	static double evaluateGTRGAMMAINVAR (int ex1, int ex2, int wptr, int iptr,
1852	double x1_start, double x2_start,
1853	double tipVector, double tFreqs, double invariants,
1854	unsigned char tipX1, int n, double diagptable, const boolean fastScaling)
1855	{
1856	int i, j, k;
1857	double x1, x2;
1858	double
1859	freqs[4],
1860	scaler = 0.25 * (1.0 - invariants),
1861	sum = 0.0,
1862	term;
1863
1864	freqs[0] = tFreqs[0] * invariants;
1865	freqs[1] = tFreqs[1] * invariants;
1866	freqs[2] = tFreqs[2] * invariants;
1867	freqs[3] = tFreqs[3] * invariants;
1868
1869	if(tipX1)
1870	{
1871	for (i = 0; i < n; i++)
1872	{
1873	x1 = &(tipVector[4 * tipX1[i]]);
1874	x2 = &x2_start[16 * i];
1875
1876	for(j = 0, term = 0.0; j < 4; j++)
1877	for(k = 0; k < 4; k++)
1878	term += x1[k] * x2[j * 4 + k] * diagptable[j * 4 + k];
1879
1880	if(iptr[i] < 4)
1881	if(fastScaling)
1882	term = LOG(((scaler * FABS(term)) + freqs[iptr[i]]));
1883	else
1884	term = LOG(((scaler * FABS(term)) + freqs[iptr[i]])) + ex2[i] * LOG(minlikelihood);
1885	else
1886	if(fastScaling)
1887	term = LOG(scaler * FABS(term));
1888	else
1889	term = LOG(scaler * FABS(term)) + (ex2[i] * LOG(minlikelihood));
1890
1891	sum += wptr[i] * term;
1892	}
1893	}
1894	else
1895	{
1896
1897	for (i = 0; i < n; i++)
1898	{
1899	x1 = &x1_start[16 * i];
1900	x2 = &x2_start[16 * i];
1901
1902	for(j = 0, term = 0.0; j < 4; j++)
1903	for(k = 0; k < 4; k++)
1904	term += x1[j * 4 + k] * x2[j * 4 + k] * diagptable[j * 4 + k];
1905
1906	if(iptr[i] < 4)
1907	if(fastScaling)
1908	term = LOG(((scaler * FABS(term)) + freqs[iptr[i]]));
1909	else
1910	term = LOG(((scaler * FABS(term)) + freqs[iptr[i]])) + (ex2[i] + ex1[i]) * LOG(minlikelihood);
1911	else
1912	if(fastScaling)
1913	term = LOG(scaler * FABS(term));
1914	else
1915	term = LOG(scaler * FABS(term)) + ((ex1[i] + ex2[i]) * LOG(minlikelihood));
1916
1917	sum += wptr[i] * term;
1918	}
1919	}
1920
1921	return sum;
1922	}
1923
1924
1925
1926
1927	static double evaluateGTRGAMMAPROT (int ex1, int ex2, int *wptr,
1928	double x1, double x2,
1929	double *tipVector,
1930	unsigned char tipX1, int n, double diagptable, const boolean fastScaling)
1931	{
1932	double sum = 0.0, term;
1933	int i, j, l;
1934	double left, right;
1935
1936	if(tipX1)
1937	{
1938	for (i = 0; i < n; i++)
1939	{
1940	#ifdef __SIM_SSE3
1941	__m128d tv = _mm_setzero_pd();
1942	left = &(tipVector[20 * tipX1[i]]);
1943
1944	for(j = 0, term = 0.0; j < 4; j++)
1945	{
1946	double d = &diagptable[j 20];
1947	right = &(x2[80 * i + 20 * j]);
1948	for(l = 0; l < 20; l+=2)
1949	{
1950	__m128d mul = _mm_mul_pd(_mm_load_pd(&left[l]), _mm_load_pd(&right[l]));
1951	tv = _mm_add_pd(tv, _mm_mul_pd(mul, _mm_load_pd(&d[l])));
1952	}
1953	}
1954	tv = _mm_hadd_pd(tv, tv);
1955	_mm_storel_pd(&term, tv);
1956
1957	#else
1958	left = &(tipVector[20 * tipX1[i]]);
1959
1960	for(j = 0, term = 0.0; j < 4; j++)
1961	{
1962	right = &(x2[80 * i + 20 * j]);
1963	for(l = 0; l < 20; l++)
1964	term += left[l] * right[l] * diagptable[j * 20 + l];
1965	}
1966	#endif
1967
1968	if(fastScaling)
1969	term = LOG(0.25 * FABS(term));
1970	else
1971	term = LOG(0.25 * FABS(term)) + (ex2[i] * LOG(minlikelihood));
1972
1973	sum += wptr[i] * term;
1974	}
1975	}
1976	else
1977	{
1978	for (i = 0; i < n; i++)
1979	{
1980	#ifdef __SIM_SSE3
1981	__m128d tv = _mm_setzero_pd();
1982
1983	for(j = 0, term = 0.0; j < 4; j++)
1984	{
1985	double d = &diagptable[j 20];
1986	left = &(x1[80 * i + 20 * j]);
1987	right = &(x2[80 * i + 20 * j]);
1988
1989	for(l = 0; l < 20; l+=2)
1990	{
1991	__m128d mul = _mm_mul_pd(_mm_load_pd(&left[l]), _mm_load_pd(&right[l]));
1992	tv = _mm_add_pd(tv, _mm_mul_pd(mul, _mm_load_pd(&d[l])));
1993	}
1994	}
1995	tv = _mm_hadd_pd(tv, tv);
1996	_mm_storel_pd(&term, tv);
1997	#else
1998	for(j = 0, term = 0.0; j < 4; j++)
1999	{
2000	left = &(x1[80 * i + 20 * j]);
2001	right = &(x2[80 * i + 20 * j]);
2002
2003	for(l = 0; l < 20; l++)
2004	term += left[l] * right[l] * diagptable[j * 20 + l];
2005	}
2006	#endif
2007
2008	if(fastScaling)
2009	term = LOG(0.25 * FABS(term));
2010	else
2011	term = LOG(0.25 * FABS(term)) + ((ex1[i] + ex2[i])*LOG(minlikelihood));
2012
2013	sum += wptr[i] * term;
2014	}
2015	}
2016
2017	return sum;
2018	}
2019
2020
2021	static double evaluateGTRGAMMAPROT_LG4(int ex1, int ex2, int *wptr,
2022	double x1, double x2,
2023	double *tipVector[4],
2024	unsigned char tipX1, int n, double diagptable, const boolean fastScaling, double *weights)
2025	{
2026	double sum = 0.0, term;
2027	int i, j, l;
2028	double left, right;
2029
2030	if(tipX1)
2031	{
2032	for (i = 0; i < n; i++)
2033	{
2034	for(j = 0, term = 0.0; j < 4; j++)
2035	{
2036	double
2037	t = 0.0;
2038
2039	left = &(tipVector[j][20 * tipX1[i]]);
2040	right = &(x2[80 * i + 20 * j]);
2041
2042	for(l = 0; l < 20; l++)
2043	t += left[l] * right[l] * diagptable[j * 20 + l];
2044
2045	term += weights[j] * t;
2046	}
2047
2048	if(fastScaling)
2049	term = LOG(FABS(term));
2050	else
2051	term = LOG(FABS(term)) + (ex2[i] * LOG(minlikelihood));
2052
2053	sum += wptr[i] * term;
2054	}
2055	}
2056	else
2057	{
2058	for (i = 0; i < n; i++)
2059	{
2060	for(j = 0, term = 0.0; j < 4; j++)
2061	{
2062	double
2063	t = 0.0;
2064
2065	left = &(x1[80 * i + 20 * j]);
2066	right = &(x2[80 * i + 20 * j]);
2067
2068	for(l = 0; l < 20; l++)
2069	t += left[l] * right[l] * diagptable[j * 20 + l];
2070
2071	term += weights[j] * t;
2072	}
2073
2074	if(fastScaling)
2075	term = LOG(FABS(term));
2076	else
2077	term = LOG(FABS(term)) + ((ex1[i] + ex2[i])*LOG(minlikelihood));
2078
2079	sum += wptr[i] * term;
2080	}
2081	}
2082
2083	return sum;
2084	}
2085
2086
2087
2088
2089	#ifdef __SIM_SSE3
2090
2091	static double evaluateGTRGAMMAPROT_GAPPED_SAVE (int ex1, int ex2, int *wptr,
2092	double x1, double x2,
2093	double *tipVector,
2094	unsigned char tipX1, int n, double diagptable, const boolean fastScaling,
2095	double x1_gapColumn, double x2_gapColumn, unsigned int x1_gap, unsigned int x2_gap)
2096	{
2097	double sum = 0.0, term;
2098	int i, j, l;
2099	double
2100	*left,
2101	*right,
2102	*x1_ptr = x1,
2103	*x2_ptr = x2,
2104	*x1v,
2105	*x2v;
2106
2107	if(tipX1)
2108	{
2109	for (i = 0; i < n; i++)
2110	{
2111	if(x2_gap[i / 32] & mask32[i % 32])
2112	x2v = x2_gapColumn;
2113	else
2114	{
2115	x2v = x2_ptr;
2116	x2_ptr += 80;
2117	}
2118
2119	__m128d tv = _mm_setzero_pd();
2120	left = &(tipVector[20 * tipX1[i]]);
2121
2122	for(j = 0, term = 0.0; j < 4; j++)
2123	{
2124	double d = &diagptable[j 20];
2125	right = &(x2v[20 * j]);
2126	for(l = 0; l < 20; l+=2)
2127	{
2128	__m128d mul = _mm_mul_pd(_mm_load_pd(&left[l]), _mm_load_pd(&right[l]));
2129	tv = _mm_add_pd(tv, _mm_mul_pd(mul, _mm_load_pd(&d[l])));
2130	}
2131	}
2132
2133	tv = _mm_hadd_pd(tv, tv);
2134	_mm_storel_pd(&term, tv);
2135
2136
2137	if(fastScaling)
2138	term = LOG(0.25 * term);
2139	else
2140	term = LOG(0.25 * term) + (ex2[i] * LOG(minlikelihood));
2141
2142	sum += wptr[i] * term;
2143	}
2144	}
2145	else
2146	{
2147	for (i = 0; i < n; i++)
2148	{
2149	if(x1_gap[i / 32] & mask32[i % 32])
2150	x1v = x1_gapColumn;
2151	else
2152	{
2153	x1v = x1_ptr;
2154	x1_ptr += 80;
2155	}
2156
2157	if(x2_gap[i / 32] & mask32[i % 32])
2158	x2v = x2_gapColumn;
2159	else
2160	{
2161	x2v = x2_ptr;
2162	x2_ptr += 80;
2163	}
2164
2165	__m128d tv = _mm_setzero_pd();
2166
2167	for(j = 0, term = 0.0; j < 4; j++)
2168	{
2169	double d = &diagptable[j 20];
2170	left = &(x1v[20 * j]);
2171	right = &(x2v[20 * j]);
2172
2173	for(l = 0; l < 20; l+=2)
2174	{
2175	__m128d mul = _mm_mul_pd(_mm_load_pd(&left[l]), _mm_load_pd(&right[l]));
2176	tv = _mm_add_pd(tv, _mm_mul_pd(mul, _mm_load_pd(&d[l])));
2177	}
2178	}
2179	tv = _mm_hadd_pd(tv, tv);
2180	_mm_storel_pd(&term, tv);
2181
2182	if(fastScaling)
2183	term = LOG(0.25 * term);
2184	else
2185	term = LOG(0.25 * term) + ((ex1[i] + ex2[i])*LOG(minlikelihood));
2186
2187	sum += wptr[i] * term;
2188	}
2189	}
2190
2191	return sum;
2192	}
2193
2194
2195	#endif
2196
2197
2198
2199
2200
2201
2202	static double evaluateGTRGAMMASECONDARY (int ex1, int ex2, int *wptr,
2203	double x1, double x2,
2204	double *tipVector,
2205	unsigned char tipX1, int n, double diagptable, const boolean fastScaling)
2206	{
2207	double sum = 0.0, term;
2208	int i, j, l;
2209	double left, right;
2210
2211	if(tipX1)
2212	{
2213	for (i = 0; i < n; i++)
2214	{
2215	left = &(tipVector[16 * tipX1[i]]);
2216
2217	for(j = 0, term = 0.0; j < 4; j++)
2218	{
2219	right = &(x2[64 * i + 16 * j]);
2220
2221	for(l = 0; l < 16; l++)
2222	term += left[l] * right[l] * diagptable[j * 16 + l];
2223	}
2224
2225	if(fastScaling)
2226	term = LOG(0.25 * FABS(term));
2227	else
2228	term = LOG(0.25 * FABS(term)) + (ex2[i] * LOG(minlikelihood));
2229
2230	sum += wptr[i] * term;
2231	}
2232	}
2233	else
2234	{
2235	for (i = 0; i < n; i++)
2236	{
2237
2238	for(j = 0, term = 0.0; j < 4; j++)
2239	{
2240	left = &(x1[64 * i + 16 * j]);
2241	right = &(x2[64 * i + 16 * j]);
2242
2243	for(l = 0; l < 16; l++)
2244	term += left[l] * right[l] * diagptable[j * 16 + l];
2245	}
2246
2247	if(fastScaling)
2248	term = LOG(0.25 * FABS(term));
2249	else
2250	term = LOG(0.25 * FABS(term)) + ((ex1[i] + ex2[i])*LOG(minlikelihood));
2251
2252	sum += wptr[i] * term;
2253	}
2254	}
2255
2256	return sum;
2257	}
2258
2259	static double evaluateGTRGAMMASECONDARY_6 (int ex1, int ex2, int *wptr,
2260	double x1, double x2,
2261	double *tipVector,
2262	unsigned char tipX1, int n, double diagptable, const boolean fastScaling)
2263	{
2264	double sum = 0.0, term;
2265	int i, j, l;
2266	double left, right;
2267
2268	if(tipX1)
2269	{
2270	for (i = 0; i < n; i++)
2271	{
2272	left = &(tipVector[6 * tipX1[i]]);
2273
2274	for(j = 0, term = 0.0; j < 4; j++)
2275	{
2276	right = &(x2[24 * i + 6 * j]);
2277
2278	for(l = 0; l < 6; l++)
2279	term += left[l] * right[l] * diagptable[j * 6 + l];
2280	}
2281
2282	if(fastScaling)
2283	term = LOG(0.25 * FABS(term));
2284	else
2285	term = LOG(0.25 * FABS(term)) + (ex2[i] * LOG(minlikelihood));
2286
2287	sum += wptr[i] * term;
2288	}
2289	}
2290	else
2291	{
2292	for (i = 0; i < n; i++)
2293	{
2294
2295	for(j = 0, term = 0.0; j < 4; j++)
2296	{
2297	left = &(x1[24 * i + 6 * j]);
2298	right = &(x2[24 * i + 6 * j]);
2299
2300	for(l = 0; l < 6; l++)
2301	term += left[l] * right[l] * diagptable[j * 6 + l];
2302	}
2303
2304	if(fastScaling)
2305	term = LOG(0.25 * FABS(term));
2306	else
2307	term = LOG(0.25 * FABS(term)) + ((ex1[i] + ex2[i])*LOG(minlikelihood));
2308
2309	sum += wptr[i] * term;
2310	}
2311	}
2312
2313	return sum;
2314	}
2315
2316	static double evaluateGTRGAMMASECONDARY_7 (int ex1, int ex2, int *wptr,
2317	double x1, double x2,
2318	double *tipVector,
2319	unsigned char tipX1, int n, double diagptable, const boolean fastScaling)
2320	{
2321	double sum = 0.0, term;
2322	int i, j, l;
2323	double left, right;
2324
2325	if(tipX1)
2326	{
2327	for (i = 0; i < n; i++)
2328	{
2329	left = &(tipVector[7 * tipX1[i]]);
2330
2331	for(j = 0, term = 0.0; j < 4; j++)
2332	{
2333	right = &(x2[28 * i + 7 * j]);
2334
2335	for(l = 0; l < 7; l++)
2336	term += left[l] * right[l] * diagptable[j * 7 + l];
2337	}
2338
2339	if(fastScaling)
2340	term = LOG(0.25 * FABS(term));
2341	else
2342	term = LOG(0.25 * FABS(term)) + (ex2[i] * LOG(minlikelihood));
2343
2344	sum += wptr[i] * term;
2345	}
2346	}
2347	else
2348	{
2349	for (i = 0; i < n; i++)
2350	{
2351
2352	for(j = 0, term = 0.0; j < 4; j++)
2353	{
2354	left = &(x1[28 * i + 7 * j]);
2355	right = &(x2[28 * i + 7 * j]);
2356
2357	for(l = 0; l < 7; l++)
2358	term += left[l] * right[l] * diagptable[j * 7 + l];
2359	}
2360
2361	if(fastScaling)
2362	term = LOG(0.25 * FABS(term));
2363	else
2364	term = LOG(0.25 * FABS(term)) + ((ex1[i] + ex2[i])*LOG(minlikelihood));
2365
2366	sum += wptr[i] * term;
2367	}
2368	}
2369
2370	return sum;
2371	}
2372
2373	static double evaluateGTRGAMMAPROTINVAR (int ex1, int ex2, int wptr, int iptr,
2374	double x1, double x2,
2375	double tipVector,double tFreqs, double invariants,
2376	unsigned char tipX1, int n, double diagptable, const boolean fastScaling)
2377	{
2378	double
2379	sum = 0.0, term, freqs[20],
2380	scaler = 0.25 * (1.0 - invariants);
2381	int i, j, l;
2382	double left, right;
2383
2384	for(i = 0; i < 20; i++)
2385	freqs[i] = tFreqs[i] * invariants;
2386
2387	if(tipX1)
2388	{
2389	for (i = 0; i < n; i++)
2390	{
2391	left = &(tipVector[20 * tipX1[i]]);
2392
2393	for(j = 0, term = 0.0; j < 4; j++)
2394	{
2395	right = &(x2[80 * i + 20 * j]);
2396
2397	for(l = 0; l < 20; l++)
2398	term += left[l] * right[l] * diagptable[j * 20 + l];
2399	}
2400
2401	if(iptr[i] < 20)
2402	if(fastScaling)
2403	term = LOG(((scaler * FABS(term)) + freqs[iptr[i]]));
2404	else
2405	term = LOG(((scaler * FABS(term)) + freqs[iptr[i]])) + ex2[i] * LOG(minlikelihood);
2406	else
2407	if(fastScaling)
2408	term = LOG(scaler * FABS(term));
2409	else
2410	term = LOG(scaler * FABS(term)) + (ex2[i] * LOG(minlikelihood));
2411
2412	sum += wptr[i] * term;
2413	}
2414	}
2415	else
2416	{
2417	for (i = 0; i < n; i++)
2418	{
2419	for(j = 0, term = 0.0; j < 4; j++)
2420	{
2421	left = &(x1[80 * i + 20 * j]);
2422	right = &(x2[80 * i + 20 * j]);
2423
2424	for(l = 0; l < 20; l++)
2425	term += left[l] * right[l] * diagptable[j * 20 + l];
2426	}
2427
2428	if(iptr[i] < 20)
2429	if(fastScaling)
2430	term = LOG(((scaler * FABS(term)) + freqs[iptr[i]]));
2431	else
2432	term = LOG(((scaler * FABS(term)) + freqs[iptr[i]])) + (ex1[i] + ex2[i]) * LOG(minlikelihood);
2433	else
2434	if(fastScaling)
2435	term = LOG(scaler * FABS(term));
2436	else
2437	term = LOG(scaler * FABS(term)) + ((ex1[i] + ex2[i]) * LOG(minlikelihood));
2438	sum += wptr[i] * term;
2439	}
2440	}
2441
2442	return sum;
2443	}
2444
2445	static double evaluateGTRGAMMASECONDARYINVAR (int ex1, int ex2, int wptr, int iptr,
2446	double x1, double x2,
2447	double tipVector,double tFreqs, double invariants,
2448	unsigned char tipX1, int n, double diagptable, const boolean fastScaling)
2449	{
2450	double
2451	sum = 0.0, term, freqs[16],
2452	scaler = 0.25 * (1.0 - invariants);
2453	int i, j, l;
2454	double left, right;
2455
2456	for(i = 0; i < 16; i++)
2457	freqs[i] = tFreqs[i] * invariants;
2458
2459	if(tipX1)
2460	{
2461	for (i = 0; i < n; i++)
2462	{
2463	left = &(tipVector[16 * tipX1[i]]);
2464
2465	for(j = 0, term = 0.0; j < 4; j++)
2466	{
2467	right = &(x2[64 * i + 16 * j]);
2468
2469	for(l = 0; l < 16; l++)
2470	term += left[l] * right[l] * diagptable[j * 16 + l];
2471	}
2472
2473	if(iptr[i] < 16)
2474	if(fastScaling)
2475	term = LOG(((scaler * FABS(term)) + freqs[iptr[i]]));
2476	else
2477	term = LOG(((scaler * FABS(term)) + freqs[iptr[i]])) + ex2[i] * LOG(minlikelihood);
2478	else
2479	if(fastScaling)
2480	term = LOG(scaler * FABS(term));
2481	else
2482	term = LOG(scaler * FABS(term)) + (ex2[i] * LOG(minlikelihood));
2483
2484	sum += wptr[i] * term;
2485	}
2486	}
2487	else
2488	{
2489	for (i = 0; i < n; i++)
2490	{
2491	for(j = 0, term = 0.0; j < 4; j++)
2492	{
2493	left = &(x1[64 * i + 16 * j]);
2494	right = &(x2[64 * i + 16 * j]);
2495
2496	for(l = 0; l < 16; l++)
2497	term += left[l] * right[l] * diagptable[j * 16 + l];
2498	}
2499
2500	if(iptr[i] < 16)
2501	if(fastScaling)
2502	term = LOG(((scaler * FABS(term)) + freqs[iptr[i]]));
2503	else
2504	term = LOG(((scaler * FABS(term)) + freqs[iptr[i]])) + (ex1[i] + ex2[i]) * LOG(minlikelihood);
2505	else
2506	if(fastScaling)
2507	term = LOG(scaler * FABS(term));
2508	else
2509	term = LOG(scaler * FABS(term)) + (ex1[i] + ex2[i]) * LOG(minlikelihood);
2510
2511	sum += wptr[i] * term;
2512	}
2513	}
2514
2515	return sum;
2516	}
2517
2518	static double evaluateGTRGAMMASECONDARYINVAR_6 (int ex1, int ex2, int wptr, int iptr,
2519	double x1, double x2,
2520	double tipVector,double tFreqs, double invariants,
2521	unsigned char tipX1, int n, double diagptable, const boolean fastScaling)
2522	{
2523	double
2524	sum = 0.0, term, freqs[6],
2525	scaler = 0.25 * (1.0 - invariants);
2526	int i, j, l;
2527	double left, right;
2528
2529	for(i = 0; i < 6; i++)
2530	freqs[i] = tFreqs[i] * invariants;
2531
2532	if(tipX1)
2533	{
2534	for (i = 0; i < n; i++)
2535	{
2536	left = &(tipVector[6 * tipX1[i]]);
2537
2538	for(j = 0, term = 0.0; j < 4; j++)
2539	{
2540	right = &(x2[24 * i + 6 * j]);
2541
2542	for(l = 0; l < 6; l++)
2543	term += left[l] * right[l] * diagptable[j * 6 + l];
2544	}
2545
2546	if(iptr[i] < 6)
2547	if(fastScaling)
2548	term = LOG(((scaler * FABS(term)) + freqs[iptr[i]]));
2549	else
2550	term = LOG(((scaler * FABS(term)) + freqs[iptr[i]])) + ex2[i] * LOG(minlikelihood);
2551	else
2552	if(fastScaling)
2553	term = LOG(scaler * FABS(term));
2554	else
2555	term = LOG(scaler * FABS(term)) + (ex2[i] * LOG(minlikelihood));
2556
2557	sum += wptr[i] * term;
2558	}
2559	}
2560	else
2561	{
2562	for (i = 0; i < n; i++)
2563	{
2564	for(j = 0, term = 0.0; j < 4; j++)
2565	{
2566	left = &(x1[24 * i + 6 * j]);
2567	right = &(x2[24 * i + 6 * j]);
2568
2569	for(l = 0; l < 6; l++)
2570	term += left[l] * right[l] * diagptable[j * 6 + l];
2571	}
2572
2573	if(iptr[i] < 6)
2574	if(fastScaling)
2575	term = LOG(((scaler * FABS(term)) + freqs[iptr[i]]));
2576	else
2577	term = LOG(((scaler * FABS(term)) + freqs[iptr[i]])) + (ex2[i] + ex1[i]) * LOG(minlikelihood);
2578	else
2579	if(fastScaling)
2580	term = LOG(scaler * FABS(term));
2581	else
2582	term = LOG(scaler * FABS(term)) + ((ex1[i] + ex2[i]) * LOG(minlikelihood));
2583
2584	sum += wptr[i] * term;
2585	}
2586	}
2587
2588	return sum;
2589	}
2590
2591	static double evaluateGTRGAMMASECONDARYINVAR_7 (int ex1, int ex2, int wptr, int iptr,
2592	double x1, double x2,
2593	double tipVector,double tFreqs, double invariants,
2594	unsigned char tipX1, int n, double diagptable, const boolean fastScaling)
2595	{
2596	double
2597	sum = 0.0, term, freqs[7],
2598	scaler = 0.25 * (1.0 - invariants);
2599	int i, j, l;
2600	double left, right;
2601
2602	for(i = 0; i < 7; i++)
2603	freqs[i] = tFreqs[i] * invariants;
2604
2605	if(tipX1)
2606	{
2607	for (i = 0; i < n; i++)
2608	{
2609	left = &(tipVector[7 * tipX1[i]]);
2610
2611	for(j = 0, term = 0.0; j < 4; j++)
2612	{
2613	right = &(x2[28 * i + 7 * j]);
2614
2615	for(l = 0; l < 7; l++)
2616	term += left[l] * right[l] * diagptable[j * 7 + l];
2617	}
2618
2619	if(iptr[i] < 7)
2620	if(fastScaling)
2621	term = LOG(((scaler * FABS(term)) + freqs[iptr[i]]));
2622	else
2623	term = LOG(((scaler * FABS(term)) + freqs[iptr[i]])) + ex2[i] * LOG(minlikelihood);
2624	else
2625	if(fastScaling)
2626	term = LOG(scaler * FABS(term));
2627	else
2628	term = LOG(scaler * FABS(term)) + (ex2[i] * LOG(minlikelihood));
2629
2630	sum += wptr[i] * term;
2631	}
2632	}
2633	else
2634	{
2635	for (i = 0; i < n; i++)
2636	{
2637	for(j = 0, term = 0.0; j < 4; j++)
2638	{
2639	left = &(x1[28 * i + 7 * j]);
2640	right = &(x2[28 * i + 7 * j]);
2641
2642	for(l = 0; l < 7; l++)
2643	term += left[l] * right[l] * diagptable[j * 7 + l];
2644	}
2645
2646	if(iptr[i] < 7)
2647	if(fastScaling)
2648	term = LOG(((scaler * FABS(term)) + freqs[iptr[i]]));
2649	else
2650	term = LOG(((scaler * FABS(term)) + freqs[iptr[i]])) + (ex2[i] + ex1[i]) * LOG(minlikelihood);
2651	else
2652	if(fastScaling)
2653	term = LOG(scaler * FABS(term));
2654	else
2655	term = LOG(scaler * FABS(term)) + ((ex1[i] + ex2[i]) * LOG(minlikelihood));
2656
2657	sum += wptr[i] * term;
2658	}
2659	}
2660
2661	return sum;
2662	}
2663
2664
2665	double evaluateIterative(tree *tr, boolean writeVector)
2666	{
2667	double
2668	*pz = tr->td[0].ti[0].qz,
2669	result = 0.0;
2670
2671	int
2672	rateHet,
2673	pNumber = tr->td[0].ti[0].pNumber,
2674	qNumber = tr->td[0].ti[0].qNumber,
2675	model;
2676
2677	if(tr->rateHetModel == CAT)
2678	rateHet = 1;
2679	else
2680	rateHet = 4;
2681
2682	newviewIterative(tr);
2683
2684	if(writeVector)
2685	assert(!tr->useFastScaling);
2686
2687	#ifdef _DEBUG_MULTI_EPA
2688	printf("EV: ");
2689	#endif
2690
2691	for(model = 0; model < tr->NumberOfModels; model++)
2692	{
2693	#ifdef _DEBUG_MULTI_EPA
2694	printf("%d ", tr->executeModel[model]);
2695	#endif
2696
2697	if(tr->executeModel[model])
2698	{
2699	int
2700	width = tr->partitionData[model].width,
2701	states = tr->partitionData[model].states;
2702
2703	double
2704	z,
2705	partitionLikelihood = 0.0,
2706	*_vector;
2707
2708	int
2709	ex1 = (int)NULL,
2710	ex2 = (int)NULL;
2711
2712	unsigned int
2713	x1_gap = (unsigned int)NULL,
2714	x2_gap = (unsigned int)NULL;
2715
2716	double
2717	*weights = tr->partitionData[model].weights,
2718	x1_start = (double)NULL,
2719	x2_start = (double)NULL,
2720	diagptable = (double)NULL,
2721	x1_gapColumn = (double)NULL,
2722	x2_gapColumn = (double)NULL;
2723
2724	unsigned char
2725	tip = (unsigned char)NULL;
2726
2727	if(writeVector)
2728	_vector = tr->partitionData[model].perSiteLL;
2729	else
2730	_vector = (double*)NULL;
2731
2732
2733	diagptable = tr->partitionData[model].left;
2734
2735
2736	if(isTip(pNumber, tr->mxtips) \|\| isTip(qNumber, tr->mxtips))
2737	{
2738	if(isTip(qNumber, tr->mxtips))
2739	{
2740	x2_start = tr->partitionData[model].xVector[pNumber - tr->mxtips -1];
2741
2742	if(!tr->useFastScaling)
2743	ex2 = tr->partitionData[model].expVector[pNumber - tr->mxtips - 1];
2744
2745	if(tr->saveMemory)
2746	{
2747	x2_gap = &(tr->partitionData[model].gapVector[pNumber * tr->partitionData[model].gapVectorLength]);
2748	x2_gapColumn = &tr->partitionData[model].gapColumn[(pNumber - tr->mxtips - 1) * states * rateHet];
2749	}
2750
2751	tip = tr->partitionData[model].yVector[qNumber];
2752	}
2753	else
2754	{
2755
2756	x2_start = tr->partitionData[model].xVector[qNumber - tr->mxtips - 1];
2757
2758
2759	if(!tr->useFastScaling)
2760	ex2 = tr->partitionData[model].expVector[qNumber - tr->mxtips - 1];
2761
2762	if(tr->saveMemory)
2763	{
2764	x2_gap = &(tr->partitionData[model].gapVector[qNumber * tr->partitionData[model].gapVectorLength]);
2765	x2_gapColumn = &tr->partitionData[model].gapColumn[(qNumber - tr->mxtips - 1) * states * rateHet];
2766	}
2767
2768	tip = tr->partitionData[model].yVector[pNumber];
2769	}
2770	}
2771	else
2772	{
2773	if(tr->saveMemory)
2774	{
2775	x1_gap = &(tr->partitionData[model].gapVector[pNumber * tr->partitionData[model].gapVectorLength]);
2776	x2_gap = &(tr->partitionData[model].gapVector[qNumber * tr->partitionData[model].gapVectorLength]);
2777	x1_gapColumn = &tr->partitionData[model].gapColumn[(pNumber - tr->mxtips - 1) * states * rateHet];
2778	x2_gapColumn = &tr->partitionData[model].gapColumn[(qNumber - tr->mxtips - 1) * states * rateHet];
2779	}
2780
2781	x1_start = tr->partitionData[model].xVector[pNumber - tr->mxtips - 1];
2782	x2_start = tr->partitionData[model].xVector[qNumber - tr->mxtips - 1];
2783
2784	if(!tr->useFastScaling)
2785	{
2786	ex1 = tr->partitionData[model].expVector[pNumber - tr->mxtips - 1];
2787	ex2 = tr->partitionData[model].expVector[qNumber - tr->mxtips - 1];
2788	}
2789	}
2790
2791
2792	if(tr->multiBranch)
2793	z = pz[model];
2794	else
2795	z = pz[0];
2796
2797	if(writeVector)
2798	{
2799	switch(tr->rateHetModel)
2800	{
2801	case CAT:
2802	{
2803	calcDiagptableFlex(z, tr->partitionData[model].numberOfCategories, tr->partitionData[model].perSiteRates, tr->partitionData[model].EIGN, diagptable, states);
2804
2805	partitionLikelihood = evaluateCatFlex(ex1, ex2, tr->partitionData[model].rateCategory, tr->partitionData[model].wgt,
2806	x1_start, x2_start, tr->partitionData[model].tipVector,
2807	tip, width, diagptable, _vector, writeVector, tr->useFastScaling, states);
2808	}
2809	break;
2810	case GAMMA:
2811	{
2812	if(tr->partitionData[model].protModels == LG4 \|\| tr->partitionData[model].protModels == LG4X)
2813	{
2814	calcDiagptableFlex_LG4(z, 4, tr->partitionData[model].gammaRates, tr->partitionData[model].EIGN_LG4, diagptable, 20);
2815
2816
2817	partitionLikelihood = evaluateGammaFlex_LG4(ex1, ex2, tr->partitionData[model].wgt,
2818	x1_start, x2_start, tr->partitionData[model].tipVector_LG4,
2819	tip, width, diagptable, _vector, writeVector, tr->useFastScaling, states, weights);
2820	}
2821	else
2822	{
2823	calcDiagptableFlex(z, 4, tr->partitionData[model].gammaRates, tr->partitionData[model].EIGN, diagptable, states);
2824
2825	partitionLikelihood = evaluateGammaFlex(ex1, ex2, tr->partitionData[model].wgt,
2826	x1_start, x2_start, tr->partitionData[model].tipVector,
2827	tip, width, diagptable, _vector, writeVector, tr->useFastScaling, states);
2828	}
2829	}
2830	break;
2831	case GAMMA_I:
2832	{
2833	calcDiagptableFlex(z, 4, tr->partitionData[model].gammaRates, tr->partitionData[model].EIGN, diagptable, states);
2834
2835	partitionLikelihood = evaluateGammaInvarFlex(ex1, ex2, tr->partitionData[model].wgt, tr->partitionData[model].invariant,
2836	x1_start, x2_start,
2837	tr->partitionData[model].tipVector, tr->partitionData[model].frequencies,
2838	tr->partitionData[model].propInvariant,
2839	tip, width, diagptable, _vector, writeVector, tr->useFastScaling, states);
2840	}
2841	break;
2842	default:
2843	assert(0);
2844	}
2845	}
2846	else
2847	{
2848	switch(tr->partitionData[model].dataType)
2849	{
2850	case BINARY_DATA:
2851	switch(tr->rateHetModel)
2852	{
2853	case CAT:
2854	{
2855	calcDiagptable(z, BINARY_DATA, tr->partitionData[model].numberOfCategories, tr->partitionData[model].perSiteRates, tr->partitionData[model].EIGN, diagptable);
2856
2857	partitionLikelihood = evaluateGTRCAT_BINARY(ex1, ex2, tr->partitionData[model].rateCategory, tr->partitionData[model].wgt,
2858	x1_start, x2_start, tr->partitionData[model].tipVector,
2859	tip, width, diagptable, tr->useFastScaling);
2860	}
2861	break;
2862	case GAMMA:
2863	{
2864	calcDiagptable(z, BINARY_DATA, 4, tr->partitionData[model].gammaRates, tr->partitionData[model].EIGN, diagptable);
2865
2866	partitionLikelihood = evaluateGTRGAMMA_BINARY(ex1, ex2, tr->partitionData[model].wgt,
2867	x1_start, x2_start, tr->partitionData[model].tipVector,
2868	tip, width, diagptable, tr->useFastScaling);
2869	}
2870	break;
2871	case GAMMA_I:
2872	{
2873	calcDiagptable(z, BINARY_DATA, 4, tr->partitionData[model].gammaRates, tr->partitionData[model].EIGN, diagptable);
2874
2875	partitionLikelihood = evaluateGTRGAMMAINVAR_BINARY(ex1, ex2, tr->partitionData[model].wgt, tr->partitionData[model].invariant,
2876	x1_start, x2_start,
2877	tr->partitionData[model].tipVector, tr->partitionData[model].frequencies,
2878	tr->partitionData[model].propInvariant,
2879	tip, width, diagptable, tr->useFastScaling);
2880	}
2881	break;
2882	default:
2883	assert(0);
2884	}
2885	break;
2886	case DNA_DATA:
2887	switch(tr->rateHetModel)
2888	{
2889	case CAT:
2890
2891	calcDiagptable(z, DNA_DATA, tr->partitionData[model].numberOfCategories, tr->partitionData[model].perSiteRates, tr->partitionData[model].EIGN, diagptable);
2892	#ifdef __SIM_SSE3
2893	if(tr->saveMemory)
2894	{
2895	partitionLikelihood = evaluateGTRCAT_SAVE(ex1, ex2, tr->partitionData[model].rateCategory, tr->partitionData[model].wgt,
2896	x1_start, x2_start, tr->partitionData[model].tipVector,
2897	tip, width, diagptable, tr->useFastScaling, x1_gapColumn, x2_gapColumn, x1_gap, x2_gap);
2898	}
2899	else
2900	#endif
2901	partitionLikelihood = evaluateGTRCAT(ex1, ex2, tr->partitionData[model].rateCategory, tr->partitionData[model].wgt,
2902	x1_start, x2_start, tr->partitionData[model].tipVector,
2903	tip, width, diagptable, tr->useFastScaling);
2904	break;
2905	case GAMMA:
2906
2907	calcDiagptable(z, DNA_DATA, 4, tr->partitionData[model].gammaRates, tr->partitionData[model].EIGN, diagptable);
2908	#ifdef __SIM_SSE3
2909	if(tr->saveMemory)
2910	partitionLikelihood = evaluateGTRGAMMA_GAPPED_SAVE(ex1, ex2, tr->partitionData[model].wgt,
2911	x1_start, x2_start, tr->partitionData[model].tipVector,
2912	tip, width, diagptable, tr->useFastScaling,
2913	x1_gapColumn, x2_gapColumn, x1_gap, x2_gap);
2914	else
2915	#endif
2916
2917	partitionLikelihood = evaluateGTRGAMMA(ex1, ex2, tr->partitionData[model].wgt,
2918	x1_start, x2_start, tr->partitionData[model].tipVector,
2919	tip, width, diagptable, tr->useFastScaling);
2920	break;
2921	case GAMMA_I:
2922	{
2923	calcDiagptable(z, DNA_DATA, 4, tr->partitionData[model].gammaRates, tr->partitionData[model].EIGN, diagptable);
2924
2925	partitionLikelihood = evaluateGTRGAMMAINVAR(ex1, ex2, tr->partitionData[model].wgt, tr->partitionData[model].invariant,
2926	x1_start, x2_start,
2927	tr->partitionData[model].tipVector, tr->partitionData[model].frequencies,
2928	tr->partitionData[model].propInvariant,
2929	tip, width, diagptable, tr->useFastScaling);
2930	}
2931	break;
2932	default:
2933	assert(0);
2934	}
2935	break;
2936	case AA_DATA:
2937	switch(tr->rateHetModel)
2938	{
2939	case CAT:
2940
2941	calcDiagptable(z, AA_DATA, tr->partitionData[model].numberOfCategories, tr->partitionData[model].perSiteRates, tr->partitionData[model].EIGN, diagptable);
2942	#ifdef __SIM_SSE3
2943	if(tr->saveMemory)
2944	{
2945	partitionLikelihood = evaluateGTRCATPROT_SAVE(ex1, ex2, tr->partitionData[model].rateCategory, tr->partitionData[model].wgt,
2946	x1_start, x2_start, tr->partitionData[model].tipVector,
2947	tip, width, diagptable, tr->useFastScaling, x1_gapColumn, x2_gapColumn, x1_gap, x2_gap);
2948	}
2949	else
2950	#endif
2951	partitionLikelihood = evaluateGTRCATPROT(ex1, ex2, tr->partitionData[model].rateCategory, tr->partitionData[model].wgt,
2952	x1_start, x2_start, tr->partitionData[model].tipVector,
2953	tip, width, diagptable, tr->useFastScaling);
2954
2955	break;
2956	case GAMMA:
2957	if(tr->partitionData[model].protModels == LG4 \|\| tr->partitionData[model].protModels == LG4X)
2958	{
2959	calcDiagptableFlex_LG4(z, 4, tr->partitionData[model].gammaRates, tr->partitionData[model].EIGN_LG4, diagptable, 20);
2960
2961	partitionLikelihood = evaluateGTRGAMMAPROT_LG4(ex1, ex2, tr->partitionData[model].wgt,
2962	x1_start, x2_start, tr->partitionData[model].tipVector_LG4,
2963	tip, width, diagptable, tr->useFastScaling, weights);
2964
2965	}
2966	else
2967	{
2968	calcDiagptable(z, AA_DATA, 4, tr->partitionData[model].gammaRates, tr->partitionData[model].EIGN, diagptable);
2969	#ifdef __SIM_SSE3
2970	if(tr->saveMemory)
2971	partitionLikelihood = evaluateGTRGAMMAPROT_GAPPED_SAVE(ex1, ex2, tr->partitionData[model].wgt,
2972	x1_start, x2_start, tr->partitionData[model].tipVector,
2973	tip, width, diagptable, tr->useFastScaling,
2974	x1_gapColumn, x2_gapColumn, x1_gap, x2_gap);
2975	else
2976	#endif
2977	partitionLikelihood = evaluateGTRGAMMAPROT(ex1, ex2, tr->partitionData[model].wgt,
2978	x1_start, x2_start, tr->partitionData[model].tipVector,
2979	tip, width, diagptable, tr->useFastScaling);
2980	}
2981	break;
2982	case GAMMA_I:
2983	{
2984	calcDiagptable(z, AA_DATA, 4, tr->partitionData[model].gammaRates, tr->partitionData[model].EIGN, diagptable);
2985
2986	partitionLikelihood = evaluateGTRGAMMAPROTINVAR(ex1, ex2, tr->partitionData[model].wgt, tr->partitionData[model].invariant,
2987	x1_start, x2_start,
2988	tr->partitionData[model].tipVector, tr->partitionData[model].frequencies,
2989	tr->partitionData[model].propInvariant,
2990	tip, width, diagptable, tr->useFastScaling);
2991	}
2992	break;
2993	default:
2994	assert(0);
2995	}
2996	break;
2997	case GENERIC_32:
2998	switch(tr->rateHetModel)
2999	{
3000	case CAT:
3001	{
3002	calcDiagptableFlex(z, tr->partitionData[model].numberOfCategories, tr->partitionData[model].perSiteRates, tr->partitionData[model].EIGN, diagptable, states);
3003
3004	partitionLikelihood = evaluateCatFlex(ex1, ex2, tr->partitionData[model].rateCategory, tr->partitionData[model].wgt,
3005	x1_start, x2_start, tr->partitionData[model].tipVector,
3006	tip, width, diagptable, _vector, writeVector, tr->useFastScaling, states);
3007	}
3008	break;
3009	case GAMMA:
3010	{
3011	calcDiagptableFlex(z, 4, tr->partitionData[model].gammaRates, tr->partitionData[model].EIGN, diagptable, states);
3012
3013	partitionLikelihood = evaluateGammaFlex(ex1, ex2, tr->partitionData[model].wgt,
3014	x1_start, x2_start, tr->partitionData[model].tipVector,
3015	tip, width, diagptable, _vector, writeVector, tr->useFastScaling, states);
3016	}
3017	break;
3018	case GAMMA_I:
3019	{
3020	calcDiagptableFlex(z, 4, tr->partitionData[model].gammaRates, tr->partitionData[model].EIGN, diagptable, states);
3021
3022	partitionLikelihood = evaluateGammaInvarFlex(ex1, ex2, tr->partitionData[model].wgt, tr->partitionData[model].invariant,
3023	x1_start, x2_start,
3024	tr->partitionData[model].tipVector, tr->partitionData[model].frequencies,
3025	tr->partitionData[model].propInvariant,
3026	tip, width, diagptable, _vector, writeVector, tr->useFastScaling, states);
3027	}
3028	break;
3029	default:
3030	assert(0);
3031	}
3032	break;
3033	case SECONDARY_DATA:
3034	switch(tr->rateHetModel)
3035	{
3036	case CAT:
3037	{
3038	calcDiagptable(z, SECONDARY_DATA, tr->partitionData[model].numberOfCategories, tr->partitionData[model].perSiteRates, tr->partitionData[model].EIGN, diagptable);
3039
3040	partitionLikelihood = evaluateGTRCATSECONDARY(ex1, ex2, tr->partitionData[model].rateCategory, tr->partitionData[model].wgt,
3041	x1_start, x2_start, tr->partitionData[model].tipVector,
3042	tip, width, diagptable, tr->useFastScaling);
3043	}
3044	break;
3045	case GAMMA:
3046	{
3047	calcDiagptable(z, SECONDARY_DATA, 4, tr->partitionData[model].gammaRates, tr->partitionData[model].EIGN, diagptable);
3048
3049	partitionLikelihood = evaluateGTRGAMMASECONDARY(ex1, ex2, tr->partitionData[model].wgt,
3050	x1_start, x2_start, tr->partitionData[model].tipVector,
3051	tip, width, diagptable, tr->useFastScaling);
3052	}
3053	break;
3054	case GAMMA_I:
3055	{
3056	calcDiagptable(z, SECONDARY_DATA, 4, tr->partitionData[model].gammaRates, tr->partitionData[model].EIGN, diagptable);
3057
3058	partitionLikelihood = evaluateGTRGAMMASECONDARYINVAR(ex1, ex2, tr->partitionData[model].wgt, tr->partitionData[model].invariant,
3059	x1_start, x2_start,
3060	tr->partitionData[model].tipVector, tr->partitionData[model].frequencies,
3061	tr->partitionData[model].propInvariant,
3062	tip, width, diagptable, tr->useFastScaling);
3063	}
3064	break;
3065	default:
3066	assert(0);
3067	}
3068	break;
3069	case SECONDARY_DATA_6:
3070	switch(tr->rateHetModel)
3071	{
3072	case CAT:
3073	{
3074	calcDiagptable(z, SECONDARY_DATA_6, tr->partitionData[model].numberOfCategories, tr->partitionData[model].perSiteRates, tr->partitionData[model].EIGN, diagptable);
3075
3076	partitionLikelihood = evaluateGTRCATSECONDARY_6(ex1, ex2, tr->partitionData[model].rateCategory, tr->partitionData[model].wgt,
3077	x1_start, x2_start, tr->partitionData[model].tipVector,
3078	tip, width, diagptable, tr->useFastScaling);
3079	}
3080	break;
3081	case GAMMA:
3082	{
3083	calcDiagptable(z, SECONDARY_DATA_6, 4, tr->partitionData[model].gammaRates, tr->partitionData[model].EIGN, diagptable);
3084
3085	partitionLikelihood = evaluateGTRGAMMASECONDARY_6(ex1, ex2, tr->partitionData[model].wgt,
3086	x1_start, x2_start, tr->partitionData[model].tipVector,
3087	tip, width, diagptable, tr->useFastScaling);
3088	}
3089	break;
3090	case GAMMA_I:
3091	{
3092	calcDiagptable(z, SECONDARY_DATA_6, 4, tr->partitionData[model].gammaRates, tr->partitionData[model].EIGN, diagptable);
3093
3094	partitionLikelihood = evaluateGTRGAMMASECONDARYINVAR_6(ex1, ex2, tr->partitionData[model].wgt, tr->partitionData[model].invariant,
3095	x1_start, x2_start,
3096	tr->partitionData[model].tipVector, tr->partitionData[model].frequencies,
3097	tr->partitionData[model].propInvariant,
3098	tip, width, diagptable, tr->useFastScaling);
3099	}
3100	break;
3101	default:
3102	assert(0);
3103	}
3104	break;
3105	case SECONDARY_DATA_7:
3106	switch(tr->rateHetModel)
3107	{
3108	case CAT:
3109	{
3110	calcDiagptable(z, SECONDARY_DATA_7, tr->partitionData[model].numberOfCategories, tr->partitionData[model].perSiteRates, tr->partitionData[model].EIGN, diagptable);
3111
3112	partitionLikelihood = evaluateGTRCATSECONDARY_7(ex1, ex2, tr->partitionData[model].rateCategory, tr->partitionData[model].wgt,
3113	x1_start, x2_start, tr->partitionData[model].tipVector,
3114	tip, width, diagptable, tr->useFastScaling);
3115	}
3116	break;
3117	case GAMMA:
3118	{
3119	calcDiagptable(z, SECONDARY_DATA_7, 4, tr->partitionData[model].gammaRates, tr->partitionData[model].EIGN, diagptable);
3120
3121	partitionLikelihood = evaluateGTRGAMMASECONDARY_7(ex1, ex2, tr->partitionData[model].wgt,
3122	x1_start, x2_start, tr->partitionData[model].tipVector,
3123	tip, width, diagptable, tr->useFastScaling);
3124	}
3125	break;
3126	case GAMMA_I:
3127	{
3128	calcDiagptable(z, SECONDARY_DATA_7, 4, tr->partitionData[model].gammaRates, tr->partitionData[model].EIGN, diagptable);
3129
3130	partitionLikelihood = evaluateGTRGAMMASECONDARYINVAR_7(ex1, ex2, tr->partitionData[model].wgt, tr->partitionData[model].invariant,
3131	x1_start, x2_start,
3132	tr->partitionData[model].tipVector, tr->partitionData[model].frequencies,
3133	tr->partitionData[model].propInvariant,
3134	tip, width, diagptable, tr->useFastScaling);
3135	}
3136	break;
3137	default:
3138	assert(0);
3139	}
3140	break;
3141	default:
3142	assert(0);
3143	}
3144	}
3145
3146	if(width > 0)
3147	{
3148	assert(partitionLikelihood < 0.0);
3149
3150	if(tr->useFastScaling)
3151	partitionLikelihood += (tr->partitionData[model].globalScaler[pNumber] + tr->partitionData[model].globalScaler[qNumber]) * LOG(minlikelihood);
3152	}
3153
3154	result += partitionLikelihood;
3155	tr->perPartitionLH[model] = partitionLikelihood;
3156	}
3157	}
3158	#ifdef _DEBUG_MULTI_EPA
3159	printf("\n");
3160	#endif
3161	return result;
3162	}
3163
3164
3165
3166
3167	double evaluateGeneric (tree *tr, nodeptr p)
3168	{
3169	volatile
3170	double result;
3171
3172	nodeptr
3173	q = p->back;
3174
3175	int
3176	i;
3177
3178
3179	tr->td[0].ti[0].pNumber = p->number;
3180	tr->td[0].ti[0].qNumber = q->number;
3181
3182	for(i = 0; i < tr->numBranches; i++)
3183	tr->td[0].ti[0].qz[i] = q->z[i];
3184
3185	tr->td[0].count = 1;
3186
3187	if(!p->x)
3188	computeTraversalInfo(p, &(tr->td[0].ti[0]), &(tr->td[0].count), tr->mxtips, tr->numBranches);
3189
3190	if(!q->x)
3191	computeTraversalInfo(q, &(tr->td[0].ti[0]), &(tr->td[0].count), tr->mxtips, tr->numBranches);
3192
3193	#ifdef _USE_PTHREADS
3194	{
3195	int j;
3196
3197	masterBarrier(THREAD_EVALUATE, tr);
3198
3199	if(tr->NumberOfModels == 1)
3200	{
3201	for(i = 0, result = 0.0; i < NumberOfThreads; i++)
3202	result += reductionBuffer[i];
3203
3204	tr->perPartitionLH[0] = result;
3205	}
3206	else
3207	{
3208	volatile
3209	double partitionResult;
3210
3211	result = 0.0;
3212
3213	for(j = 0; j < tr->NumberOfModels; j++)
3214	{
3215	for(i = 0, partitionResult = 0.0; i < NumberOfThreads; i++)
3216	partitionResult += reductionBuffer[i * tr->NumberOfModels + j];
3217	result += partitionResult;
3218	tr->perPartitionLH[j] = partitionResult;
3219	}
3220	}
3221	}
3222	#else
3223	result = evaluateIterative(tr, FALSE);
3224	#endif
3225
3226	tr->likelihood = result;
3227
3228	return result;
3229	}
3230
3231
3232
3233
3234	double evaluateGenericInitrav (tree *tr, nodeptr p)
3235	{
3236	volatile double
3237	result;
3238
3239	determineFullTraversal(p, tr);
3240
3241	#ifdef _USE_PTHREADS
3242	{
3243	int
3244	i,
3245	j;
3246
3247	masterBarrier(THREAD_EVALUATE, tr);
3248
3249	if(tr->NumberOfModels == 1)
3250	{
3251	for(i = 0, result = 0.0; i < NumberOfThreads; i++)
3252	result += reductionBuffer[i];
3253
3254	tr->perPartitionLH[0] = result;
3255	}
3256	else
3257	{
3258	volatile double
3259	partitionResult;
3260
3261	result = 0.0;
3262
3263	for(j = 0; j < tr->NumberOfModels; j++)
3264	{
3265	for(i = 0, partitionResult = 0.0; i < NumberOfThreads; i++)
3266	partitionResult += reductionBuffer[i * tr->NumberOfModels + j];
3267	result += partitionResult;
3268	tr->perPartitionLH[j] = partitionResult;
3269	}
3270	}
3271	}
3272	#else
3273	result = evaluateIterative(tr, FALSE);
3274	#endif
3275
3276	tr->likelihood = result;
3277
3278	return result;
3279	}
3280
3281
3282
3283
3284	void onlyInitrav(tree *tr, nodeptr p)
3285	{
3286	determineFullTraversal(p, tr);
3287
3288	#ifdef _USE_PTHREADS
3289	masterBarrier(THREAD_NEWVIEW, tr);
3290	#else
3291	newviewIterative(tr);
3292	#endif
3293	}
3294
3295
3296
3297
3298
3299
3300	#ifdef _USE_PTHREADS
3301
3302	double evalCL(tree tr, double x2, int _ex2, unsigned char _tip, double *pz, int insertion)
3303	{
3304	double
3305	x1_start = (double)NULL,
3306	result = 0.0;
3307
3308	int
3309	ex1 = (int)NULL,
3310	model,
3311	columnCounter,
3312	offsetCounter;
3313
3314	unsigned char
3315	tip = (unsigned char)NULL;
3316
3317	setPartitionMask(tr, insertion, tr->executeModel);
3318
3319	#ifdef _DEBUG_MULTI_EPA
3320	if(tr->threadID == THREAD_TO_DEBUG)
3321	printf("EV %s: ", tr->nameList[tr->inserts[insertion]]);
3322	#endif
3323
3324	for(model = 0, columnCounter = 0, offsetCounter = 0; model < tr->NumberOfModels; model++)
3325	{
3326	int
3327	width = tr->partitionData[model].upper - tr->partitionData[model].lower;
3328
3329	#ifdef _DEBUG_MULTI_EPA
3330	if(tr->threadID == THREAD_TO_DEBUG)
3331	printf("%d", tr->executeModel[model]);
3332	#endif
3333
3334	if(tr->executeModel[model])
3335	{
3336	int
3337	*ex2,
3338	*rateCategory,
3339	*wgt,
3340	*invariant;
3341
3342	double
3343	*x2_start,
3344	z,
3345	partitionLikelihood,
3346	*diagptable = tr->partitionData[model].left;
3347
3348
3349	rateCategory = &tr->contiguousRateCategory[columnCounter];
3350	wgt = &tr->contiguousWgt[columnCounter];
3351	invariant = &tr->contiguousInvariant[columnCounter];
3352	tip = &_tip[columnCounter];
3353	x2_start = &x2[offsetCounter];
3354	ex2 = &_ex2[columnCounter];
3355
3356	if(tr->multiBranch)
3357	z = pz[model];
3358	else
3359	z = pz[0];
3360
3361	switch(tr->partitionData[model].dataType)
3362	{
3363	case BINARY_DATA:
3364	switch(tr->rateHetModel)
3365	{
3366	case CAT:
3367	calcDiagptable(z, BINARY_DATA, tr->partitionData[model].numberOfCategories, tr->partitionData[model].perSiteRates, tr->partitionData[model].EIGN, diagptable);
3368
3369	partitionLikelihood = evaluateGTRCAT_BINARY(ex1, ex2, rateCategory, wgt,
3370	x1_start, x2_start, tr->partitionData[model].tipVector,
3371	tip, width, diagptable, tr->useFastScaling);
3372	break;
3373	case GAMMA:
3374	calcDiagptable(z, BINARY_DATA, 4, tr->partitionData[model].gammaRates, tr->partitionData[model].EIGN, diagptable);
3375
3376	partitionLikelihood = evaluateGTRGAMMA_BINARY(ex1, ex2,wgt,
3377	x1_start, x2_start, tr->partitionData[model].tipVector,
3378	tip, width, diagptable, tr->useFastScaling);
3379
3380	break;
3381	case GAMMA_I:
3382	calcDiagptable(z, BINARY_DATA, 4, tr->partitionData[model].gammaRates, tr->partitionData[model].EIGN, diagptable);
3383
3384	partitionLikelihood = evaluateGTRGAMMAINVAR_BINARY(ex1, ex2,wgt, invariant,
3385	x1_start, x2_start,
3386	tr->partitionData[model].tipVector, tr->partitionData[model].frequencies,
3387	tr->partitionData[model].propInvariant,
3388	tip, width, diagptable, tr->useFastScaling);
3389	break;
3390	default:
3391	assert(0);
3392	}
3393	break;
3394	case DNA_DATA:
3395	switch(tr->rateHetModel)
3396	{
3397	case CAT:
3398	calcDiagptable(z, DNA_DATA, tr->partitionData[model].numberOfCategories, tr->partitionData[model].perSiteRates, tr->partitionData[model].EIGN, diagptable);
3399
3400	partitionLikelihood = evaluateGTRCAT(ex1, ex2, rateCategory,wgt,
3401	x1_start, x2_start, tr->partitionData[model].tipVector,
3402	tip, width, diagptable, tr->useFastScaling);
3403	break;
3404	case GAMMA:
3405	calcDiagptable(z, DNA_DATA, 4, tr->partitionData[model].gammaRates, tr->partitionData[model].EIGN, diagptable);
3406
3407	partitionLikelihood = evaluateGTRGAMMA(ex1, ex2,wgt,
3408	x1_start, x2_start, tr->partitionData[model].tipVector,
3409	tip, width, diagptable, tr->useFastScaling);
3410	break;
3411	case GAMMA_I:
3412	calcDiagptable(z, DNA_DATA, 4, tr->partitionData[model].gammaRates, tr->partitionData[model].EIGN, diagptable);
3413
3414	partitionLikelihood = evaluateGTRGAMMAINVAR(ex1, ex2,wgt,invariant,
3415	x1_start, x2_start,
3416	tr->partitionData[model].tipVector, tr->partitionData[model].frequencies,
3417	tr->partitionData[model].propInvariant,
3418	tip, width, diagptable, tr->useFastScaling);
3419	break;
3420	default:
3421	assert(0);
3422	}
3423	break;
3424	case AA_DATA:
3425	switch(tr->rateHetModel)
3426	{
3427	case CAT:
3428	calcDiagptable(z, AA_DATA, tr->partitionData[model].numberOfCategories, tr->partitionData[model].perSiteRates, tr->partitionData[model].EIGN, diagptable);
3429
3430	partitionLikelihood = evaluateGTRCATPROT(ex1, ex2, rateCategory,wgt,
3431	x1_start, x2_start, tr->partitionData[model].tipVector,
3432	tip, width, diagptable, tr->useFastScaling);
3433	break;
3434	case GAMMA:
3435	calcDiagptable(z, AA_DATA, 4, tr->partitionData[model].gammaRates, tr->partitionData[model].EIGN, diagptable);
3436
3437	partitionLikelihood = evaluateGTRGAMMAPROT(ex1, ex2,wgt,
3438	x1_start, x2_start, tr->partitionData[model].tipVector,
3439	tip, width, diagptable, tr->useFastScaling);
3440	break;
3441	case GAMMA_I:
3442	calcDiagptable(z, AA_DATA, 4, tr->partitionData[model].gammaRates, tr->partitionData[model].EIGN, diagptable);
3443
3444	partitionLikelihood = evaluateGTRGAMMAPROTINVAR(ex1, ex2,wgt,invariant,
3445	x1_start, x2_start,
3446	tr->partitionData[model].tipVector, tr->partitionData[model].frequencies,
3447	tr->partitionData[model].propInvariant,
3448	tip, width, diagptable, tr->useFastScaling);
3449	break;
3450	default:
3451	assert(0);
3452	}
3453	break;
3454	case SECONDARY_DATA:
3455	switch(tr->rateHetModel)
3456	{
3457	case CAT:
3458	calcDiagptable(z, SECONDARY_DATA, tr->partitionData[model].numberOfCategories, tr->partitionData[model].perSiteRates, tr->partitionData[model].EIGN, diagptable);
3459
3460	partitionLikelihood = evaluateGTRCATSECONDARY(ex1, ex2, rateCategory,wgt,
3461	x1_start, x2_start, tr->partitionData[model].tipVector,
3462	tip, width, diagptable, tr->useFastScaling);
3463	break;
3464	case GAMMA:
3465	calcDiagptable(z, SECONDARY_DATA, 4, tr->partitionData[model].gammaRates, tr->partitionData[model].EIGN, diagptable);
3466
3467	partitionLikelihood = evaluateGTRGAMMASECONDARY(ex1, ex2,wgt,
3468	x1_start, x2_start, tr->partitionData[model].tipVector,
3469	tip, width, diagptable, tr->useFastScaling);
3470	break;
3471	case GAMMA_I:
3472	calcDiagptable(z, SECONDARY_DATA, 4, tr->partitionData[model].gammaRates, tr->partitionData[model].EIGN, diagptable);
3473
3474	partitionLikelihood = evaluateGTRGAMMASECONDARYINVAR(ex1, ex2,wgt,invariant,
3475	x1_start, x2_start,
3476	tr->partitionData[model].tipVector, tr->partitionData[model].frequencies,
3477	tr->partitionData[model].propInvariant,
3478	tip, width, diagptable, tr->useFastScaling);
3479	break;
3480	default:
3481	assert(0);
3482	}
3483	break;
3484	case SECONDARY_DATA_6:
3485	switch(tr->rateHetModel)
3486	{
3487	case CAT:
3488	calcDiagptable(z, SECONDARY_DATA_6, tr->partitionData[model].numberOfCategories, tr->partitionData[model].perSiteRates, tr->partitionData[model].EIGN, diagptable);
3489
3490	partitionLikelihood = evaluateGTRCATSECONDARY_6(ex1, ex2, rateCategory,wgt,
3491	x1_start, x2_start, tr->partitionData[model].tipVector,
3492	tip, width, diagptable, tr->useFastScaling);
3493	break;
3494	case GAMMA:
3495	calcDiagptable(z, SECONDARY_DATA_6, 4, tr->partitionData[model].gammaRates, tr->partitionData[model].EIGN, diagptable);
3496
3497	partitionLikelihood = evaluateGTRGAMMASECONDARY_6(ex1, ex2,wgt,
3498	x1_start, x2_start, tr->partitionData[model].tipVector,
3499	tip, width, diagptable, tr->useFastScaling);
3500	break;
3501	case GAMMA_I:
3502	calcDiagptable(z, SECONDARY_DATA_6, 4, tr->partitionData[model].gammaRates, tr->partitionData[model].EIGN, diagptable);
3503
3504	partitionLikelihood = evaluateGTRGAMMASECONDARYINVAR_6(ex1, ex2,wgt,invariant,
3505	x1_start, x2_start,
3506	tr->partitionData[model].tipVector, tr->partitionData[model].frequencies,
3507	tr->partitionData[model].propInvariant,
3508	tip, width, diagptable, tr->useFastScaling);
3509	break;
3510	default:
3511	assert(0);
3512	}
3513	break;
3514	case SECONDARY_DATA_7:
3515	switch(tr->rateHetModel)
3516	{
3517	case CAT:
3518	calcDiagptable(z, SECONDARY_DATA_7, tr->partitionData[model].numberOfCategories, tr->partitionData[model].perSiteRates, tr->partitionData[model].EIGN, diagptable);
3519
3520	partitionLikelihood = evaluateGTRCATSECONDARY_7(ex1, ex2, rateCategory,wgt,
3521	x1_start, x2_start, tr->partitionData[model].tipVector,
3522	tip, width, diagptable, tr->useFastScaling);
3523	break;
3524	case GAMMA:
3525	calcDiagptable(z, SECONDARY_DATA_7, 4, tr->partitionData[model].gammaRates, tr->partitionData[model].EIGN, diagptable);
3526
3527	partitionLikelihood = evaluateGTRGAMMASECONDARY_7(ex1, ex2,wgt,
3528	x1_start, x2_start, tr->partitionData[model].tipVector,
3529	tip, width, diagptable, tr->useFastScaling);
3530	break;
3531	case GAMMA_I:
3532	calcDiagptable(z, SECONDARY_DATA_7, 4, tr->partitionData[model].gammaRates, tr->partitionData[model].EIGN, diagptable);
3533
3534	partitionLikelihood = evaluateGTRGAMMASECONDARYINVAR_7(ex1, ex2,wgt,invariant,
3535	x1_start, x2_start,
3536	tr->partitionData[model].tipVector, tr->partitionData[model].frequencies,
3537	tr->partitionData[model].propInvariant,
3538	tip, width, diagptable, tr->useFastScaling);
3539	break;
3540	default:
3541	assert(0);
3542	}
3543	break;
3544	default:
3545	assert(0);
3546	}
3547
3548	assert(!tr->useFastScaling);
3549
3550	/* error ? */
3551
3552	tr->perPartitionLH[model] = partitionLikelihood;
3553
3554	result += partitionLikelihood;
3555	}
3556
3557	columnCounter += width;
3558	offsetCounter += width * tr->partitionData[model].states * tr->discreteRateCategories;
3559	}
3560
3561	resetPartitionMask(tr, tr->executeModel);
3562	#ifdef _DEBUG_MULTI_EPA
3563	if(tr->threadID == THREAD_TO_DEBUG)
3564	printf("\n");
3565	#endif
3566	if(tr->perPartitionEPA)
3567	{
3568
3569	return (tr->perPartitionLH[tr->readPartition[insertion]]);
3570	}
3571	else
3572	{
3573	return result;
3574	}
3575	}
3576
3577
3578
3579
3580
3581
3582
3583	#endif
3584
3585
3586
3587
3588	/*****************************************************************************************************/
3589
3590
3591
3592	double evaluateGenericVector (tree *tr, nodeptr p)
3593	{
3594	volatile double result;
3595	nodeptr q = p->back;
3596	int i;
3597
3598
3599	{
3600	tr->td[0].ti[0].pNumber = p->number;
3601	tr->td[0].ti[0].qNumber = q->number;
3602
3603	for(i = 0; i < tr->numBranches; i++)
3604	tr->td[0].ti[0].qz[i] = q->z[i];
3605
3606	tr->td[0].count = 1;
3607	if(!p->x)
3608	computeTraversalInfo(p, &(tr->td[0].ti[0]), &(tr->td[0].count), tr->mxtips, tr->numBranches);
3609	if(!q->x)
3610	computeTraversalInfo(q, &(tr->td[0].ti[0]), &(tr->td[0].count), tr->mxtips, tr->numBranches);
3611
3612	#ifdef _USE_PTHREADS
3613	{
3614	int j;
3615
3616	masterBarrier(THREAD_EVALUATE_VECTOR, tr);
3617	if(tr->NumberOfModels == 1)
3618	{
3619	for(i = 0, result = 0.0; i < NumberOfThreads; i++)
3620	result += reductionBuffer[i];
3621
3622	tr->perPartitionLH[0] = result;
3623	}
3624	else
3625	{
3626	volatile double partitionResult;
3627
3628	result = 0.0;
3629
3630	for(j = 0; j < tr->NumberOfModels; j++)
3631	{
3632	for(i = 0, partitionResult = 0.0; i < NumberOfThreads; i++)
3633	partitionResult += reductionBuffer[i * tr->NumberOfModels + j];
3634	result += partitionResult;
3635	tr->perPartitionLH[j] = partitionResult;
3636	}
3637	}
3638	}
3639	#else
3640	result = evaluateIterative(tr, TRUE);
3641	#endif
3642	}
3643
3644	tr->likelihood = result;
3645
3646	return result;
3647	}

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

Download in other formats: