Context Navigation

← Previous Revision
Next Revision →
Blame
Revision Log

evaluateGenericSpecial.c

Visit:

Last change on this file was 19480, checked in by westram, 2 years ago
reintegrates 'gde' into 'trunk' adds: log:branches/gde@19460:19479
File size: 95.4 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	#if defined(__SIM_SSE3)
2672	int
2673	rateHet;
2674	#endif
2675
2676	int
2677	pNumber = tr->td[0].ti[0].pNumber,
2678	qNumber = tr->td[0].ti[0].qNumber,
2679	model;
2680
2681	#if defined(__SIM_SSE3)
2682	if(tr->rateHetModel == CAT)
2683	rateHet = 1;
2684	else
2685	rateHet = 4;
2686	#endif
2687
2688	newviewIterative(tr);
2689
2690	if(writeVector)
2691	assert(!tr->useFastScaling);
2692
2693	#ifdef _DEBUG_MULTI_EPA
2694	printf("EV: ");
2695	#endif
2696
2697	for(model = 0; model < tr->NumberOfModels; model++)
2698	{
2699	#ifdef _DEBUG_MULTI_EPA
2700	printf("%d ", tr->executeModel[model]);
2701	#endif
2702
2703	if(tr->executeModel[model])
2704	{
2705	int
2706	width = tr->partitionData[model].width,
2707	states = tr->partitionData[model].states;
2708
2709	double
2710	z,
2711	partitionLikelihood = 0.0,
2712	*_vector;
2713
2714	int
2715	ex1 = (int)NULL,
2716	ex2 = (int)NULL;
2717
2718	#if defined(__SIM_SSE3)
2719	unsigned int
2720	x1_gap = (unsigned int)NULL,
2721	x2_gap = (unsigned int)NULL;
2722	#endif
2723
2724	double
2725	*weights = tr->partitionData[model].weights,
2726	x1_start = (double)NULL,
2727	x2_start = (double)NULL,
2728	diagptable = (double)NULL;
2729
2730	#if defined(__SIM_SSE3)
2731	double
2732	x1_gapColumn = (double)NULL,
2733	x2_gapColumn = (double)NULL;
2734	#endif
2735
2736	unsigned char
2737	tip = (unsigned char)NULL;
2738
2739	if(writeVector)
2740	_vector = tr->partitionData[model].perSiteLL;
2741	else
2742	_vector = (double*)NULL;
2743
2744
2745	diagptable = tr->partitionData[model].left;
2746
2747
2748	if(isTip(pNumber, tr->mxtips) \|\| isTip(qNumber, tr->mxtips))
2749	{
2750	if(isTip(qNumber, tr->mxtips))
2751	{
2752	x2_start = tr->partitionData[model].xVector[pNumber - tr->mxtips -1];
2753
2754	if(!tr->useFastScaling)
2755	ex2 = tr->partitionData[model].expVector[pNumber - tr->mxtips - 1];
2756
2757	#if defined(__SIM_SSE3)
2758	if(tr->saveMemory)
2759	{
2760	x2_gap = &(tr->partitionData[model].gapVector[pNumber * tr->partitionData[model].gapVectorLength]);
2761	x2_gapColumn = &tr->partitionData[model].gapColumn[(pNumber - tr->mxtips - 1) * states * rateHet];
2762	}
2763	#endif
2764
2765	tip = tr->partitionData[model].yVector[qNumber];
2766	}
2767	else
2768	{
2769
2770	x2_start = tr->partitionData[model].xVector[qNumber - tr->mxtips - 1];
2771
2772
2773	if(!tr->useFastScaling)
2774	ex2 = tr->partitionData[model].expVector[qNumber - tr->mxtips - 1];
2775
2776	#if defined(__SIM_SSE3)
2777	if(tr->saveMemory)
2778	{
2779	x2_gap = &(tr->partitionData[model].gapVector[qNumber * tr->partitionData[model].gapVectorLength]);
2780	x2_gapColumn = &tr->partitionData[model].gapColumn[(qNumber - tr->mxtips - 1) * states * rateHet];
2781	}
2782	#endif
2783
2784	tip = tr->partitionData[model].yVector[pNumber];
2785	}
2786	}
2787	else
2788	{
2789	#if defined(__SIM_SSE3)
2790	if(tr->saveMemory)
2791	{
2792	x1_gap = &(tr->partitionData[model].gapVector[pNumber * tr->partitionData[model].gapVectorLength]);
2793	x2_gap = &(tr->partitionData[model].gapVector[qNumber * tr->partitionData[model].gapVectorLength]);
2794	x1_gapColumn = &tr->partitionData[model].gapColumn[(pNumber - tr->mxtips - 1) * states * rateHet];
2795	x2_gapColumn = &tr->partitionData[model].gapColumn[(qNumber - tr->mxtips - 1) * states * rateHet];
2796	}
2797	#endif
2798
2799	x1_start = tr->partitionData[model].xVector[pNumber - tr->mxtips - 1];
2800	x2_start = tr->partitionData[model].xVector[qNumber - tr->mxtips - 1];
2801
2802	if(!tr->useFastScaling)
2803	{
2804	ex1 = tr->partitionData[model].expVector[pNumber - tr->mxtips - 1];
2805	ex2 = tr->partitionData[model].expVector[qNumber - tr->mxtips - 1];
2806	}
2807	}
2808
2809
2810	if(tr->multiBranch)
2811	z = pz[model];
2812	else
2813	z = pz[0];
2814
2815	if(writeVector)
2816	{
2817	switch(tr->rateHetModel)
2818	{
2819	case CAT:
2820	{
2821	calcDiagptableFlex(z, tr->partitionData[model].numberOfCategories, tr->partitionData[model].perSiteRates, tr->partitionData[model].EIGN, diagptable, states);
2822
2823	partitionLikelihood = evaluateCatFlex(ex1, ex2, tr->partitionData[model].rateCategory, tr->partitionData[model].wgt,
2824	x1_start, x2_start, tr->partitionData[model].tipVector,
2825	tip, width, diagptable, _vector, writeVector, tr->useFastScaling, states);
2826	}
2827	break;
2828	case GAMMA:
2829	{
2830	if(tr->partitionData[model].protModels == LG4 \|\| tr->partitionData[model].protModels == LG4X)
2831	{
2832	calcDiagptableFlex_LG4(z, 4, tr->partitionData[model].gammaRates, tr->partitionData[model].EIGN_LG4, diagptable, 20);
2833
2834
2835	partitionLikelihood = evaluateGammaFlex_LG4(ex1, ex2, tr->partitionData[model].wgt,
2836	x1_start, x2_start, tr->partitionData[model].tipVector_LG4,
2837	tip, width, diagptable, _vector, writeVector, tr->useFastScaling, states, weights);
2838	}
2839	else
2840	{
2841	calcDiagptableFlex(z, 4, tr->partitionData[model].gammaRates, tr->partitionData[model].EIGN, diagptable, states);
2842
2843	partitionLikelihood = evaluateGammaFlex(ex1, ex2, tr->partitionData[model].wgt,
2844	x1_start, x2_start, tr->partitionData[model].tipVector,
2845	tip, width, diagptable, _vector, writeVector, tr->useFastScaling, states);
2846	}
2847	}
2848	break;
2849	case GAMMA_I:
2850	{
2851	calcDiagptableFlex(z, 4, tr->partitionData[model].gammaRates, tr->partitionData[model].EIGN, diagptable, states);
2852
2853	partitionLikelihood = evaluateGammaInvarFlex(ex1, ex2, tr->partitionData[model].wgt, tr->partitionData[model].invariant,
2854	x1_start, x2_start,
2855	tr->partitionData[model].tipVector, tr->partitionData[model].frequencies,
2856	tr->partitionData[model].propInvariant,
2857	tip, width, diagptable, _vector, writeVector, tr->useFastScaling, states);
2858	}
2859	break;
2860	default:
2861	assert(0);
2862	}
2863	}
2864	else
2865	{
2866	switch(tr->partitionData[model].dataType)
2867	{
2868	case BINARY_DATA:
2869	switch(tr->rateHetModel)
2870	{
2871	case CAT:
2872	{
2873	calcDiagptable(z, BINARY_DATA, tr->partitionData[model].numberOfCategories, tr->partitionData[model].perSiteRates, tr->partitionData[model].EIGN, diagptable);
2874
2875	partitionLikelihood = evaluateGTRCAT_BINARY(ex1, ex2, tr->partitionData[model].rateCategory, tr->partitionData[model].wgt,
2876	x1_start, x2_start, tr->partitionData[model].tipVector,
2877	tip, width, diagptable, tr->useFastScaling);
2878	}
2879	break;
2880	case GAMMA:
2881	{
2882	calcDiagptable(z, BINARY_DATA, 4, tr->partitionData[model].gammaRates, tr->partitionData[model].EIGN, diagptable);
2883
2884	partitionLikelihood = evaluateGTRGAMMA_BINARY(ex1, ex2, tr->partitionData[model].wgt,
2885	x1_start, x2_start, tr->partitionData[model].tipVector,
2886	tip, width, diagptable, tr->useFastScaling);
2887	}
2888	break;
2889	case GAMMA_I:
2890	{
2891	calcDiagptable(z, BINARY_DATA, 4, tr->partitionData[model].gammaRates, tr->partitionData[model].EIGN, diagptable);
2892
2893	partitionLikelihood = evaluateGTRGAMMAINVAR_BINARY(ex1, ex2, tr->partitionData[model].wgt, tr->partitionData[model].invariant,
2894	x1_start, x2_start,
2895	tr->partitionData[model].tipVector, tr->partitionData[model].frequencies,
2896	tr->partitionData[model].propInvariant,
2897	tip, width, diagptable, tr->useFastScaling);
2898	}
2899	break;
2900	default:
2901	assert(0);
2902	}
2903	break;
2904	case DNA_DATA:
2905	switch(tr->rateHetModel)
2906	{
2907	case CAT:
2908
2909	calcDiagptable(z, DNA_DATA, tr->partitionData[model].numberOfCategories, tr->partitionData[model].perSiteRates, tr->partitionData[model].EIGN, diagptable);
2910	#ifdef __SIM_SSE3
2911	if(tr->saveMemory)
2912	{
2913	partitionLikelihood = evaluateGTRCAT_SAVE(ex1, ex2, tr->partitionData[model].rateCategory, tr->partitionData[model].wgt,
2914	x1_start, x2_start, tr->partitionData[model].tipVector,
2915	tip, width, diagptable, tr->useFastScaling, x1_gapColumn, x2_gapColumn, x1_gap, x2_gap);
2916	}
2917	else
2918	#endif
2919	partitionLikelihood = evaluateGTRCAT(ex1, ex2, tr->partitionData[model].rateCategory, tr->partitionData[model].wgt,
2920	x1_start, x2_start, tr->partitionData[model].tipVector,
2921	tip, width, diagptable, tr->useFastScaling);
2922	break;
2923	case GAMMA:
2924
2925	calcDiagptable(z, DNA_DATA, 4, tr->partitionData[model].gammaRates, tr->partitionData[model].EIGN, diagptable);
2926	#ifdef __SIM_SSE3
2927	if(tr->saveMemory)
2928	partitionLikelihood = evaluateGTRGAMMA_GAPPED_SAVE(ex1, ex2, tr->partitionData[model].wgt,
2929	x1_start, x2_start, tr->partitionData[model].tipVector,
2930	tip, width, diagptable, tr->useFastScaling,
2931	x1_gapColumn, x2_gapColumn, x1_gap, x2_gap);
2932	else
2933	#endif
2934
2935	partitionLikelihood = evaluateGTRGAMMA(ex1, ex2, tr->partitionData[model].wgt,
2936	x1_start, x2_start, tr->partitionData[model].tipVector,
2937	tip, width, diagptable, tr->useFastScaling);
2938	break;
2939	case GAMMA_I:
2940	{
2941	calcDiagptable(z, DNA_DATA, 4, tr->partitionData[model].gammaRates, tr->partitionData[model].EIGN, diagptable);
2942
2943	partitionLikelihood = evaluateGTRGAMMAINVAR(ex1, ex2, tr->partitionData[model].wgt, tr->partitionData[model].invariant,
2944	x1_start, x2_start,
2945	tr->partitionData[model].tipVector, tr->partitionData[model].frequencies,
2946	tr->partitionData[model].propInvariant,
2947	tip, width, diagptable, tr->useFastScaling);
2948	}
2949	break;
2950	default:
2951	assert(0);
2952	}
2953	break;
2954	case AA_DATA:
2955	switch(tr->rateHetModel)
2956	{
2957	case CAT:
2958
2959	calcDiagptable(z, AA_DATA, tr->partitionData[model].numberOfCategories, tr->partitionData[model].perSiteRates, tr->partitionData[model].EIGN, diagptable);
2960	#ifdef __SIM_SSE3
2961	if(tr->saveMemory)
2962	{
2963	partitionLikelihood = evaluateGTRCATPROT_SAVE(ex1, ex2, tr->partitionData[model].rateCategory, tr->partitionData[model].wgt,
2964	x1_start, x2_start, tr->partitionData[model].tipVector,
2965	tip, width, diagptable, tr->useFastScaling, x1_gapColumn, x2_gapColumn, x1_gap, x2_gap);
2966	}
2967	else
2968	#endif
2969	partitionLikelihood = evaluateGTRCATPROT(ex1, ex2, tr->partitionData[model].rateCategory, tr->partitionData[model].wgt,
2970	x1_start, x2_start, tr->partitionData[model].tipVector,
2971	tip, width, diagptable, tr->useFastScaling);
2972
2973	break;
2974	case GAMMA:
2975	if(tr->partitionData[model].protModels == LG4 \|\| tr->partitionData[model].protModels == LG4X)
2976	{
2977	calcDiagptableFlex_LG4(z, 4, tr->partitionData[model].gammaRates, tr->partitionData[model].EIGN_LG4, diagptable, 20);
2978
2979	partitionLikelihood = evaluateGTRGAMMAPROT_LG4(ex1, ex2, tr->partitionData[model].wgt,
2980	x1_start, x2_start, tr->partitionData[model].tipVector_LG4,
2981	tip, width, diagptable, tr->useFastScaling, weights);
2982
2983	}
2984	else
2985	{
2986	calcDiagptable(z, AA_DATA, 4, tr->partitionData[model].gammaRates, tr->partitionData[model].EIGN, diagptable);
2987	#ifdef __SIM_SSE3
2988	if(tr->saveMemory)
2989	partitionLikelihood = evaluateGTRGAMMAPROT_GAPPED_SAVE(ex1, ex2, tr->partitionData[model].wgt,
2990	x1_start, x2_start, tr->partitionData[model].tipVector,
2991	tip, width, diagptable, tr->useFastScaling,
2992	x1_gapColumn, x2_gapColumn, x1_gap, x2_gap);
2993	else
2994	#endif
2995	partitionLikelihood = evaluateGTRGAMMAPROT(ex1, ex2, tr->partitionData[model].wgt,
2996	x1_start, x2_start, tr->partitionData[model].tipVector,
2997	tip, width, diagptable, tr->useFastScaling);
2998	}
2999	break;
3000	case GAMMA_I:
3001	{
3002	calcDiagptable(z, AA_DATA, 4, tr->partitionData[model].gammaRates, tr->partitionData[model].EIGN, diagptable);
3003
3004	partitionLikelihood = evaluateGTRGAMMAPROTINVAR(ex1, ex2, tr->partitionData[model].wgt, tr->partitionData[model].invariant,
3005	x1_start, x2_start,
3006	tr->partitionData[model].tipVector, tr->partitionData[model].frequencies,
3007	tr->partitionData[model].propInvariant,
3008	tip, width, diagptable, tr->useFastScaling);
3009	}
3010	break;
3011	default:
3012	assert(0);
3013	}
3014	break;
3015	case GENERIC_32:
3016	switch(tr->rateHetModel)
3017	{
3018	case CAT:
3019	{
3020	calcDiagptableFlex(z, tr->partitionData[model].numberOfCategories, tr->partitionData[model].perSiteRates, tr->partitionData[model].EIGN, diagptable, states);
3021
3022	partitionLikelihood = evaluateCatFlex(ex1, ex2, tr->partitionData[model].rateCategory, tr->partitionData[model].wgt,
3023	x1_start, x2_start, tr->partitionData[model].tipVector,
3024	tip, width, diagptable, _vector, writeVector, tr->useFastScaling, states);
3025	}
3026	break;
3027	case GAMMA:
3028	{
3029	calcDiagptableFlex(z, 4, tr->partitionData[model].gammaRates, tr->partitionData[model].EIGN, diagptable, states);
3030
3031	partitionLikelihood = evaluateGammaFlex(ex1, ex2, tr->partitionData[model].wgt,
3032	x1_start, x2_start, tr->partitionData[model].tipVector,
3033	tip, width, diagptable, _vector, writeVector, tr->useFastScaling, states);
3034	}
3035	break;
3036	case GAMMA_I:
3037	{
3038	calcDiagptableFlex(z, 4, tr->partitionData[model].gammaRates, tr->partitionData[model].EIGN, diagptable, states);
3039
3040	partitionLikelihood = evaluateGammaInvarFlex(ex1, ex2, tr->partitionData[model].wgt, tr->partitionData[model].invariant,
3041	x1_start, x2_start,
3042	tr->partitionData[model].tipVector, tr->partitionData[model].frequencies,
3043	tr->partitionData[model].propInvariant,
3044	tip, width, diagptable, _vector, writeVector, tr->useFastScaling, states);
3045	}
3046	break;
3047	default:
3048	assert(0);
3049	}
3050	break;
3051	case SECONDARY_DATA:
3052	switch(tr->rateHetModel)
3053	{
3054	case CAT:
3055	{
3056	calcDiagptable(z, SECONDARY_DATA, tr->partitionData[model].numberOfCategories, tr->partitionData[model].perSiteRates, tr->partitionData[model].EIGN, diagptable);
3057
3058	partitionLikelihood = evaluateGTRCATSECONDARY(ex1, ex2, tr->partitionData[model].rateCategory, tr->partitionData[model].wgt,
3059	x1_start, x2_start, tr->partitionData[model].tipVector,
3060	tip, width, diagptable, tr->useFastScaling);
3061	}
3062	break;
3063	case GAMMA:
3064	{
3065	calcDiagptable(z, SECONDARY_DATA, 4, tr->partitionData[model].gammaRates, tr->partitionData[model].EIGN, diagptable);
3066
3067	partitionLikelihood = evaluateGTRGAMMASECONDARY(ex1, ex2, tr->partitionData[model].wgt,
3068	x1_start, x2_start, tr->partitionData[model].tipVector,
3069	tip, width, diagptable, tr->useFastScaling);
3070	}
3071	break;
3072	case GAMMA_I:
3073	{
3074	calcDiagptable(z, SECONDARY_DATA, 4, tr->partitionData[model].gammaRates, tr->partitionData[model].EIGN, diagptable);
3075
3076	partitionLikelihood = evaluateGTRGAMMASECONDARYINVAR(ex1, ex2, tr->partitionData[model].wgt, tr->partitionData[model].invariant,
3077	x1_start, x2_start,
3078	tr->partitionData[model].tipVector, tr->partitionData[model].frequencies,
3079	tr->partitionData[model].propInvariant,
3080	tip, width, diagptable, tr->useFastScaling);
3081	}
3082	break;
3083	default:
3084	assert(0);
3085	}
3086	break;
3087	case SECONDARY_DATA_6:
3088	switch(tr->rateHetModel)
3089	{
3090	case CAT:
3091	{
3092	calcDiagptable(z, SECONDARY_DATA_6, tr->partitionData[model].numberOfCategories, tr->partitionData[model].perSiteRates, tr->partitionData[model].EIGN, diagptable);
3093
3094	partitionLikelihood = evaluateGTRCATSECONDARY_6(ex1, ex2, tr->partitionData[model].rateCategory, tr->partitionData[model].wgt,
3095	x1_start, x2_start, tr->partitionData[model].tipVector,
3096	tip, width, diagptable, tr->useFastScaling);
3097	}
3098	break;
3099	case GAMMA:
3100	{
3101	calcDiagptable(z, SECONDARY_DATA_6, 4, tr->partitionData[model].gammaRates, tr->partitionData[model].EIGN, diagptable);
3102
3103	partitionLikelihood = evaluateGTRGAMMASECONDARY_6(ex1, ex2, tr->partitionData[model].wgt,
3104	x1_start, x2_start, tr->partitionData[model].tipVector,
3105	tip, width, diagptable, tr->useFastScaling);
3106	}
3107	break;
3108	case GAMMA_I:
3109	{
3110	calcDiagptable(z, SECONDARY_DATA_6, 4, tr->partitionData[model].gammaRates, tr->partitionData[model].EIGN, diagptable);
3111
3112	partitionLikelihood = evaluateGTRGAMMASECONDARYINVAR_6(ex1, ex2, tr->partitionData[model].wgt, tr->partitionData[model].invariant,
3113	x1_start, x2_start,
3114	tr->partitionData[model].tipVector, tr->partitionData[model].frequencies,
3115	tr->partitionData[model].propInvariant,
3116	tip, width, diagptable, tr->useFastScaling);
3117	}
3118	break;
3119	default:
3120	assert(0);
3121	}
3122	break;
3123	case SECONDARY_DATA_7:
3124	switch(tr->rateHetModel)
3125	{
3126	case CAT:
3127	{
3128	calcDiagptable(z, SECONDARY_DATA_7, tr->partitionData[model].numberOfCategories, tr->partitionData[model].perSiteRates, tr->partitionData[model].EIGN, diagptable);
3129
3130	partitionLikelihood = evaluateGTRCATSECONDARY_7(ex1, ex2, tr->partitionData[model].rateCategory, tr->partitionData[model].wgt,
3131	x1_start, x2_start, tr->partitionData[model].tipVector,
3132	tip, width, diagptable, tr->useFastScaling);
3133	}
3134	break;
3135	case GAMMA:
3136	{
3137	calcDiagptable(z, SECONDARY_DATA_7, 4, tr->partitionData[model].gammaRates, tr->partitionData[model].EIGN, diagptable);
3138
3139	partitionLikelihood = evaluateGTRGAMMASECONDARY_7(ex1, ex2, tr->partitionData[model].wgt,
3140	x1_start, x2_start, tr->partitionData[model].tipVector,
3141	tip, width, diagptable, tr->useFastScaling);
3142	}
3143	break;
3144	case GAMMA_I:
3145	{
3146	calcDiagptable(z, SECONDARY_DATA_7, 4, tr->partitionData[model].gammaRates, tr->partitionData[model].EIGN, diagptable);
3147
3148	partitionLikelihood = evaluateGTRGAMMASECONDARYINVAR_7(ex1, ex2, tr->partitionData[model].wgt, tr->partitionData[model].invariant,
3149	x1_start, x2_start,
3150	tr->partitionData[model].tipVector, tr->partitionData[model].frequencies,
3151	tr->partitionData[model].propInvariant,
3152	tip, width, diagptable, tr->useFastScaling);
3153	}
3154	break;
3155	default:
3156	assert(0);
3157	}
3158	break;
3159	default:
3160	assert(0);
3161	}
3162	}
3163
3164	if(width > 0)
3165	{
3166	assert(partitionLikelihood < 0.0);
3167
3168	if(tr->useFastScaling)
3169	partitionLikelihood += (tr->partitionData[model].globalScaler[pNumber] + tr->partitionData[model].globalScaler[qNumber]) * LOG(minlikelihood);
3170	}
3171
3172	result += partitionLikelihood;
3173	tr->perPartitionLH[model] = partitionLikelihood;
3174	}
3175	}
3176	#ifdef _DEBUG_MULTI_EPA
3177	printf("\n");
3178	#endif
3179	return result;
3180	}
3181
3182
3183
3184
3185	double evaluateGeneric (tree *tr, nodeptr p)
3186	{
3187	volatile
3188	double result;
3189
3190	nodeptr
3191	q = p->back;
3192
3193	int
3194	i;
3195
3196
3197	tr->td[0].ti[0].pNumber = p->number;
3198	tr->td[0].ti[0].qNumber = q->number;
3199
3200	for(i = 0; i < tr->numBranches; i++)
3201	tr->td[0].ti[0].qz[i] = q->z[i];
3202
3203	tr->td[0].count = 1;
3204
3205	if(!p->x)
3206	computeTraversalInfo(p, &(tr->td[0].ti[0]), &(tr->td[0].count), tr->mxtips, tr->numBranches);
3207
3208	if(!q->x)
3209	computeTraversalInfo(q, &(tr->td[0].ti[0]), &(tr->td[0].count), tr->mxtips, tr->numBranches);
3210
3211	#ifdef _USE_PTHREADS
3212	{
3213	int j;
3214
3215	masterBarrier(THREAD_EVALUATE, tr);
3216
3217	if(tr->NumberOfModels == 1)
3218	{
3219	for(i = 0, result = 0.0; i < NumberOfThreads; i++)
3220	result += reductionBuffer[i];
3221
3222	tr->perPartitionLH[0] = result;
3223	}
3224	else
3225	{
3226	volatile
3227	double partitionResult;
3228
3229	result = 0.0;
3230
3231	for(j = 0; j < tr->NumberOfModels; j++)
3232	{
3233	for(i = 0, partitionResult = 0.0; i < NumberOfThreads; i++)
3234	partitionResult += reductionBuffer[i * tr->NumberOfModels + j];
3235	result += partitionResult;
3236	tr->perPartitionLH[j] = partitionResult;
3237	}
3238	}
3239	}
3240	#else
3241	result = evaluateIterative(tr, FALSE);
3242	#endif
3243
3244	tr->likelihood = result;
3245
3246	return result;
3247	}
3248
3249
3250
3251
3252	double evaluateGenericInitrav (tree *tr, nodeptr p)
3253	{
3254	volatile double
3255	result;
3256
3257	determineFullTraversal(p, tr);
3258
3259	#ifdef _USE_PTHREADS
3260	{
3261	int
3262	i,
3263	j;
3264
3265	masterBarrier(THREAD_EVALUATE, tr);
3266
3267	if(tr->NumberOfModels == 1)
3268	{
3269	for(i = 0, result = 0.0; i < NumberOfThreads; i++)
3270	result += reductionBuffer[i];
3271
3272	tr->perPartitionLH[0] = result;
3273	}
3274	else
3275	{
3276	volatile double
3277	partitionResult;
3278
3279	result = 0.0;
3280
3281	for(j = 0; j < tr->NumberOfModels; j++)
3282	{
3283	for(i = 0, partitionResult = 0.0; i < NumberOfThreads; i++)
3284	partitionResult += reductionBuffer[i * tr->NumberOfModels + j];
3285	result += partitionResult;
3286	tr->perPartitionLH[j] = partitionResult;
3287	}
3288	}
3289	}
3290	#else
3291	result = evaluateIterative(tr, FALSE);
3292	#endif
3293
3294	tr->likelihood = result;
3295
3296	return result;
3297	}
3298
3299
3300
3301
3302	void onlyInitrav(tree *tr, nodeptr p)
3303	{
3304	determineFullTraversal(p, tr);
3305
3306	#ifdef _USE_PTHREADS
3307	masterBarrier(THREAD_NEWVIEW, tr);
3308	#else
3309	newviewIterative(tr);
3310	#endif
3311	}
3312
3313
3314
3315
3316
3317
3318	#ifdef _USE_PTHREADS
3319
3320	double evalCL(tree tr, double x2, int _ex2, unsigned char _tip, double *pz, int insertion)
3321	{
3322	double
3323	x1_start = (double)NULL,
3324	result = 0.0;
3325
3326	int
3327	ex1 = (int)NULL,
3328	model,
3329	columnCounter,
3330	offsetCounter;
3331
3332	unsigned char
3333	tip = (unsigned char)NULL;
3334
3335	setPartitionMask(tr, insertion, tr->executeModel);
3336
3337	#ifdef _DEBUG_MULTI_EPA
3338	if(tr->threadID == THREAD_TO_DEBUG)
3339	printf("EV %s: ", tr->nameList[tr->inserts[insertion]]);
3340	#endif
3341
3342	for(model = 0, columnCounter = 0, offsetCounter = 0; model < tr->NumberOfModels; model++)
3343	{
3344	int
3345	width = tr->partitionData[model].upper - tr->partitionData[model].lower;
3346
3347	#ifdef _DEBUG_MULTI_EPA
3348	if(tr->threadID == THREAD_TO_DEBUG)
3349	printf("%d", tr->executeModel[model]);
3350	#endif
3351
3352	if(tr->executeModel[model])
3353	{
3354	int
3355	*ex2,
3356	*rateCategory,
3357	*wgt,
3358	*invariant;
3359
3360	double
3361	*x2_start,
3362	z,
3363	partitionLikelihood,
3364	*diagptable = tr->partitionData[model].left;
3365
3366
3367	rateCategory = &tr->contiguousRateCategory[columnCounter];
3368	wgt = &tr->contiguousWgt[columnCounter];
3369	invariant = &tr->contiguousInvariant[columnCounter];
3370	tip = &_tip[columnCounter];
3371	x2_start = &x2[offsetCounter];
3372	ex2 = &_ex2[columnCounter];
3373
3374	if(tr->multiBranch)
3375	z = pz[model];
3376	else
3377	z = pz[0];
3378
3379	switch(tr->partitionData[model].dataType)
3380	{
3381	case BINARY_DATA:
3382	switch(tr->rateHetModel)
3383	{
3384	case CAT:
3385	calcDiagptable(z, BINARY_DATA, tr->partitionData[model].numberOfCategories, tr->partitionData[model].perSiteRates, tr->partitionData[model].EIGN, diagptable);
3386
3387	partitionLikelihood = evaluateGTRCAT_BINARY(ex1, ex2, rateCategory, wgt,
3388	x1_start, x2_start, tr->partitionData[model].tipVector,
3389	tip, width, diagptable, tr->useFastScaling);
3390	break;
3391	case GAMMA:
3392	calcDiagptable(z, BINARY_DATA, 4, tr->partitionData[model].gammaRates, tr->partitionData[model].EIGN, diagptable);
3393
3394	partitionLikelihood = evaluateGTRGAMMA_BINARY(ex1, ex2,wgt,
3395	x1_start, x2_start, tr->partitionData[model].tipVector,
3396	tip, width, diagptable, tr->useFastScaling);
3397
3398	break;
3399	case GAMMA_I:
3400	calcDiagptable(z, BINARY_DATA, 4, tr->partitionData[model].gammaRates, tr->partitionData[model].EIGN, diagptable);
3401
3402	partitionLikelihood = evaluateGTRGAMMAINVAR_BINARY(ex1, ex2,wgt, invariant,
3403	x1_start, x2_start,
3404	tr->partitionData[model].tipVector, tr->partitionData[model].frequencies,
3405	tr->partitionData[model].propInvariant,
3406	tip, width, diagptable, tr->useFastScaling);
3407	break;
3408	default:
3409	assert(0);
3410	}
3411	break;
3412	case DNA_DATA:
3413	switch(tr->rateHetModel)
3414	{
3415	case CAT:
3416	calcDiagptable(z, DNA_DATA, tr->partitionData[model].numberOfCategories, tr->partitionData[model].perSiteRates, tr->partitionData[model].EIGN, diagptable);
3417
3418	partitionLikelihood = evaluateGTRCAT(ex1, ex2, rateCategory,wgt,
3419	x1_start, x2_start, tr->partitionData[model].tipVector,
3420	tip, width, diagptable, tr->useFastScaling);
3421	break;
3422	case GAMMA:
3423	calcDiagptable(z, DNA_DATA, 4, tr->partitionData[model].gammaRates, tr->partitionData[model].EIGN, diagptable);
3424
3425	partitionLikelihood = evaluateGTRGAMMA(ex1, ex2,wgt,
3426	x1_start, x2_start, tr->partitionData[model].tipVector,
3427	tip, width, diagptable, tr->useFastScaling);
3428	break;
3429	case GAMMA_I:
3430	calcDiagptable(z, DNA_DATA, 4, tr->partitionData[model].gammaRates, tr->partitionData[model].EIGN, diagptable);
3431
3432	partitionLikelihood = evaluateGTRGAMMAINVAR(ex1, ex2,wgt,invariant,
3433	x1_start, x2_start,
3434	tr->partitionData[model].tipVector, tr->partitionData[model].frequencies,
3435	tr->partitionData[model].propInvariant,
3436	tip, width, diagptable, tr->useFastScaling);
3437	break;
3438	default:
3439	assert(0);
3440	}
3441	break;
3442	case AA_DATA:
3443	switch(tr->rateHetModel)
3444	{
3445	case CAT:
3446	calcDiagptable(z, AA_DATA, tr->partitionData[model].numberOfCategories, tr->partitionData[model].perSiteRates, tr->partitionData[model].EIGN, diagptable);
3447
3448	partitionLikelihood = evaluateGTRCATPROT(ex1, ex2, rateCategory,wgt,
3449	x1_start, x2_start, tr->partitionData[model].tipVector,
3450	tip, width, diagptable, tr->useFastScaling);
3451	break;
3452	case GAMMA:
3453	calcDiagptable(z, AA_DATA, 4, tr->partitionData[model].gammaRates, tr->partitionData[model].EIGN, diagptable);
3454
3455	partitionLikelihood = evaluateGTRGAMMAPROT(ex1, ex2,wgt,
3456	x1_start, x2_start, tr->partitionData[model].tipVector,
3457	tip, width, diagptable, tr->useFastScaling);
3458	break;
3459	case GAMMA_I:
3460	calcDiagptable(z, AA_DATA, 4, tr->partitionData[model].gammaRates, tr->partitionData[model].EIGN, diagptable);
3461
3462	partitionLikelihood = evaluateGTRGAMMAPROTINVAR(ex1, ex2,wgt,invariant,
3463	x1_start, x2_start,
3464	tr->partitionData[model].tipVector, tr->partitionData[model].frequencies,
3465	tr->partitionData[model].propInvariant,
3466	tip, width, diagptable, tr->useFastScaling);
3467	break;
3468	default:
3469	assert(0);
3470	}
3471	break;
3472	case SECONDARY_DATA:
3473	switch(tr->rateHetModel)
3474	{
3475	case CAT:
3476	calcDiagptable(z, SECONDARY_DATA, tr->partitionData[model].numberOfCategories, tr->partitionData[model].perSiteRates, tr->partitionData[model].EIGN, diagptable);
3477
3478	partitionLikelihood = evaluateGTRCATSECONDARY(ex1, ex2, rateCategory,wgt,
3479	x1_start, x2_start, tr->partitionData[model].tipVector,
3480	tip, width, diagptable, tr->useFastScaling);
3481	break;
3482	case GAMMA:
3483	calcDiagptable(z, SECONDARY_DATA, 4, tr->partitionData[model].gammaRates, tr->partitionData[model].EIGN, diagptable);
3484
3485	partitionLikelihood = evaluateGTRGAMMASECONDARY(ex1, ex2,wgt,
3486	x1_start, x2_start, tr->partitionData[model].tipVector,
3487	tip, width, diagptable, tr->useFastScaling);
3488	break;
3489	case GAMMA_I:
3490	calcDiagptable(z, SECONDARY_DATA, 4, tr->partitionData[model].gammaRates, tr->partitionData[model].EIGN, diagptable);
3491
3492	partitionLikelihood = evaluateGTRGAMMASECONDARYINVAR(ex1, ex2,wgt,invariant,
3493	x1_start, x2_start,
3494	tr->partitionData[model].tipVector, tr->partitionData[model].frequencies,
3495	tr->partitionData[model].propInvariant,
3496	tip, width, diagptable, tr->useFastScaling);
3497	break;
3498	default:
3499	assert(0);
3500	}
3501	break;
3502	case SECONDARY_DATA_6:
3503	switch(tr->rateHetModel)
3504	{
3505	case CAT:
3506	calcDiagptable(z, SECONDARY_DATA_6, tr->partitionData[model].numberOfCategories, tr->partitionData[model].perSiteRates, tr->partitionData[model].EIGN, diagptable);
3507
3508	partitionLikelihood = evaluateGTRCATSECONDARY_6(ex1, ex2, rateCategory,wgt,
3509	x1_start, x2_start, tr->partitionData[model].tipVector,
3510	tip, width, diagptable, tr->useFastScaling);
3511	break;
3512	case GAMMA:
3513	calcDiagptable(z, SECONDARY_DATA_6, 4, tr->partitionData[model].gammaRates, tr->partitionData[model].EIGN, diagptable);
3514
3515	partitionLikelihood = evaluateGTRGAMMASECONDARY_6(ex1, ex2,wgt,
3516	x1_start, x2_start, tr->partitionData[model].tipVector,
3517	tip, width, diagptable, tr->useFastScaling);
3518	break;
3519	case GAMMA_I:
3520	calcDiagptable(z, SECONDARY_DATA_6, 4, tr->partitionData[model].gammaRates, tr->partitionData[model].EIGN, diagptable);
3521
3522	partitionLikelihood = evaluateGTRGAMMASECONDARYINVAR_6(ex1, ex2,wgt,invariant,
3523	x1_start, x2_start,
3524	tr->partitionData[model].tipVector, tr->partitionData[model].frequencies,
3525	tr->partitionData[model].propInvariant,
3526	tip, width, diagptable, tr->useFastScaling);
3527	break;
3528	default:
3529	assert(0);
3530	}
3531	break;
3532	case SECONDARY_DATA_7:
3533	switch(tr->rateHetModel)
3534	{
3535	case CAT:
3536	calcDiagptable(z, SECONDARY_DATA_7, tr->partitionData[model].numberOfCategories, tr->partitionData[model].perSiteRates, tr->partitionData[model].EIGN, diagptable);
3537
3538	partitionLikelihood = evaluateGTRCATSECONDARY_7(ex1, ex2, rateCategory,wgt,
3539	x1_start, x2_start, tr->partitionData[model].tipVector,
3540	tip, width, diagptable, tr->useFastScaling);
3541	break;
3542	case GAMMA:
3543	calcDiagptable(z, SECONDARY_DATA_7, 4, tr->partitionData[model].gammaRates, tr->partitionData[model].EIGN, diagptable);
3544
3545	partitionLikelihood = evaluateGTRGAMMASECONDARY_7(ex1, ex2,wgt,
3546	x1_start, x2_start, tr->partitionData[model].tipVector,
3547	tip, width, diagptable, tr->useFastScaling);
3548	break;
3549	case GAMMA_I:
3550	calcDiagptable(z, SECONDARY_DATA_7, 4, tr->partitionData[model].gammaRates, tr->partitionData[model].EIGN, diagptable);
3551
3552	partitionLikelihood = evaluateGTRGAMMASECONDARYINVAR_7(ex1, ex2,wgt,invariant,
3553	x1_start, x2_start,
3554	tr->partitionData[model].tipVector, tr->partitionData[model].frequencies,
3555	tr->partitionData[model].propInvariant,
3556	tip, width, diagptable, tr->useFastScaling);
3557	break;
3558	default:
3559	assert(0);
3560	}
3561	break;
3562	default:
3563	assert(0);
3564	}
3565
3566	assert(!tr->useFastScaling);
3567
3568	/* error ? */
3569
3570	tr->perPartitionLH[model] = partitionLikelihood;
3571
3572	result += partitionLikelihood;
3573	}
3574
3575	columnCounter += width;
3576	offsetCounter += width * tr->partitionData[model].states * tr->discreteRateCategories;
3577	}
3578
3579	resetPartitionMask(tr, tr->executeModel);
3580	#ifdef _DEBUG_MULTI_EPA
3581	if(tr->threadID == THREAD_TO_DEBUG)
3582	printf("\n");
3583	#endif
3584	if(tr->perPartitionEPA)
3585	{
3586
3587	return (tr->perPartitionLH[tr->readPartition[insertion]]);
3588	}
3589	else
3590	{
3591	return result;
3592	}
3593	}
3594
3595
3596
3597
3598
3599
3600
3601	#endif
3602
3603
3604
3605
3606	/*****************************************************************************************************/
3607
3608
3609
3610	double evaluateGenericVector (tree *tr, nodeptr p)
3611	{
3612	volatile double result;
3613	nodeptr q = p->back;
3614	int i;
3615
3616
3617	{
3618	tr->td[0].ti[0].pNumber = p->number;
3619	tr->td[0].ti[0].qNumber = q->number;
3620
3621	for(i = 0; i < tr->numBranches; i++)
3622	tr->td[0].ti[0].qz[i] = q->z[i];
3623
3624	tr->td[0].count = 1;
3625	if(!p->x)
3626	computeTraversalInfo(p, &(tr->td[0].ti[0]), &(tr->td[0].count), tr->mxtips, tr->numBranches);
3627	if(!q->x)
3628	computeTraversalInfo(q, &(tr->td[0].ti[0]), &(tr->td[0].count), tr->mxtips, tr->numBranches);
3629
3630	#ifdef _USE_PTHREADS
3631	{
3632	int j;
3633
3634	masterBarrier(THREAD_EVALUATE_VECTOR, tr);
3635	if(tr->NumberOfModels == 1)
3636	{
3637	for(i = 0, result = 0.0; i < NumberOfThreads; i++)
3638	result += reductionBuffer[i];
3639
3640	tr->perPartitionLH[0] = result;
3641	}
3642	else
3643	{
3644	volatile double partitionResult;
3645
3646	result = 0.0;
3647
3648	for(j = 0; j < tr->NumberOfModels; j++)
3649	{
3650	for(i = 0, partitionResult = 0.0; i < NumberOfThreads; i++)
3651	partitionResult += reductionBuffer[i * tr->NumberOfModels + j];
3652	result += partitionResult;
3653	tr->perPartitionLH[j] = partitionResult;
3654	}
3655	}
3656	}
3657	#else
3658	result = evaluateIterative(tr, TRUE);
3659	#endif
3660	}
3661
3662	tr->likelihood = result;
3663
3664	return result;
3665	}

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

Context Navigation

source: trunk/GDE/RAxML/evaluateGenericSpecial.c

Download in other formats: