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

Last change on this file was 5262, checked in by westram, 17 years ago
  • update to RAxML 7.0.3
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 39.9 KB
Line 
1/*  RAxML-VI-HPC (version 2.2) a program for sequential and parallel estimation of phylogenetic trees
2 *  Copyright August 2006 by Alexandros Stamatakis
3 *
4 *  Partially derived from
5 *  fastDNAml, a program for estimation of phylogenetic trees from sequences by Gary J. Olsen
6 * 
7 *  and
8 *
9 *  Programs of the PHYLIP package by Joe Felsenstein.
10 
11 *  This program is free software; you may redistribute it and/or modify its
12 *  under the terms of the GNU General Public License as published by the Free
13 *  Software Foundation; either version 2 of the License, or (at your option)
14 *  any later version.
15 *
16 *  This program is distributed in the hope that it will be useful, but
17 *  WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
18 *  or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
19 *  for more details.
20 *
21 *
22 *  For any other enquiries send an Email to Alexandros Stamatakis
23 *  Alexandros.Stamatakis@epfl.ch
24 *
25 *  When publishing work that is based on the results from RAxML-VI-HPC please cite:
26 *
27 *  Alexandros Stamatakis:"RAxML-VI-HPC: maximum likelihood-based phylogenetic analyses with thousands of taxa and mixed models".
28 *  Bioinformatics 2006; doi: 10.1093/bioinformatics/btl446
29 */
30
31#ifndef WIN32
32#include <unistd.h>
33#endif
34
35#include <math.h>
36#include <time.h>
37#include <stdlib.h>
38#include <stdio.h>
39#include <ctype.h>
40#include <string.h>
41#include "axml.h"
42
43
44
45
46
47static void evaluateGTRCAT_VECTOR (int *ex2, int *cptr,
48                                   double *x2_start, double *v, double *EIGN, double *rptr, double *tipVector, 
49                                   double pz, 
50                                   char *tipX1, int lower, int n, int numberOfCategories)
51{
52  double   z, lz, ki, lz1, lz2, lz3, term;   
53  int     i; 
54  double  *diagptable, *diagptable_start;
55  double *x1, *x2;             
56 
57  z = pz;
58 
59  if (z < zmin) z = zmin;
60  lz = log(z);
61 
62  lz1 = EIGN[0] * lz;
63  lz2 = EIGN[1] * lz;
64  lz3 = EIGN[2] * lz;
65 
66  diagptable = diagptable_start = (double *)malloc(sizeof(double) * numberOfCategories * 3);   
67   
68  for(i = 0; i <  numberOfCategories; i++)
69    {   
70      ki = rptr[i];     
71      *diagptable++ = exp(ki * lz1);
72      *diagptable++ = exp(ki * lz2);
73      *diagptable++ = exp(ki * lz3);   
74    }         
75 
76
77  for (i = lower; i < n; i++) 
78    {                             
79      x1 = &(tipVector[4 * tipX1[i]]);
80      x2 = &x2_start[4 * i];
81     
82      diagptable = &diagptable_start[3 * cptr[i]];                 
83     
84      term =  x1[0] * x2[0];       
85      term += x1[1] * x2[1] * *diagptable++;     
86      term += x1[2] * x2[2] * *diagptable++;       
87      term += x1[3] * x2[3] * *diagptable++; 
88     
89      term = log(term) + (ex2[i] * log(minlikelihood));                   
90      v[i] = term;     
91    }
92  free(diagptable_start);       
93}         
94 
95static void evaluateGTRCATMULT_VECTOR (int *ex2, int *cptr, int *modelptr, 
96                                       double *x2_start, double *v, double *EIGN, double *rptr, 
97                                       double *tipVector, double *pz, 
98                                       char *tipX1, int lower, int n, int numberOfCategories, int numberOfModels, int multiBranch)
99{
100  double   z, lz = 0.0, ki, lz1, lz2, lz3, term;   
101  int     i;
102  double  *diagptable, *diagptable_start; 
103  double *x1, *x2;
104  int model, modelCounter;       
105 
106  if(!multiBranch)
107    {
108      z = pz[0]; 
109      if (z < zmin) z = zmin;
110      lz = log(z);
111    }
112   
113  diagptable = diagptable_start = (double *)malloc(sizeof(double) * numberOfCategories * 3 * numberOfModels);   
114 
115  for(modelCounter = 0; modelCounter < numberOfModels; modelCounter++)
116    {
117      if(multiBranch)
118        {
119          z = pz[modelCounter]; 
120          if (z < zmin) z = zmin;
121          lz = log(z);
122        }
123
124      lz1 = EIGN[modelCounter * 3] * lz;
125      lz2 = EIGN[modelCounter * 3 + 1] * lz;
126      lz3 = EIGN[modelCounter * 3 + 2] * lz;         
127     
128      for(i = 0; i <  numberOfCategories; i++)
129        {       
130          ki = rptr[i]; 
131          *diagptable++ = exp (ki * lz1);
132          *diagptable++ = exp (ki * lz2);
133          *diagptable++ = exp (ki * lz3);       
134        }
135    }                         
136
137   
138  for (i = lower; i < n; i++) 
139    {                   
140      model = modelptr[i];
141      x1 = &(tipVector[model * 64 + 4 * tipX1[i]]);
142      x2 = &x2_start[4 * i];
143     
144      diagptable = &diagptable_start[model * 3 *  numberOfCategories + 3 * cptr[i]];               
145     
146      term =  x1[0] * x2[0];       
147      term += x1[1] * x2[1] * *diagptable++;     
148      term += x1[2] * x2[2] * *diagptable++;       
149      term += x1[3] * x2[3] * *diagptable++; 
150     
151      term = (log(term)) + (ex2[i] * log(minlikelihood));                         
152       
153      v[i] = term;                                       
154    }
155  free(diagptable_start); 
156}
157   
158static void evaluateGTRCATPROT_VECTOR (int *ex2, int *cptr,
159                                       double *x2, double *v, double *EIGN, double *rptr, double *tipVector, double pz, 
160                                       char *tipX1, int lower, int n, int numberOfCategories)
161{
162  double   z, lz, ki, lza[19];       
163  int     i, l; 
164  double  *diagptable, *diagptable_start; 
165  double term;
166  double *left, *right;           
167   
168  z = pz;
169   
170  if (z < zmin) z = zmin;
171  lz = log(z);
172
173  for(l = 0; l < 19; l++)     
174    lza[l] = EIGN[l] * lz;     
175
176  diagptable = diagptable_start = (double *)malloc(sizeof(double) * numberOfCategories * 19);   
177 
178  for(i = 0; i <  numberOfCategories; i++)
179    {   
180      ki = rptr[i];     
181       
182      *diagptable++ = exp (ki * lza[0]);           
183      *diagptable++ = exp (ki * lza[1]);
184      *diagptable++ = exp (ki * lza[2]);           
185      *diagptable++ = exp (ki * lza[3]);
186      *diagptable++ = exp (ki * lza[4]);           
187      *diagptable++ = exp (ki * lza[5]);
188      *diagptable++ = exp (ki * lza[6]);           
189      *diagptable++ = exp (ki * lza[7]);
190      *diagptable++ = exp (ki * lza[8]);           
191      *diagptable++ = exp (ki * lza[9]);
192      *diagptable++ = exp (ki * lza[10]);         
193      *diagptable++ = exp (ki * lza[11]);
194      *diagptable++ = exp (ki * lza[12]);         
195      *diagptable++ = exp (ki * lza[13]);
196      *diagptable++ = exp (ki * lza[14]);         
197      *diagptable++ = exp (ki * lza[15]);
198      *diagptable++ = exp (ki * lza[16]);         
199      *diagptable++ = exp (ki * lza[17]);
200      *diagptable++ = exp (ki * lza[18]);                 
201    }                       
202   
203
204  for (i = lower; i < n; i++) 
205    {                             
206      left = &(tipVector[20 * tipX1[i]]);
207      right = &(x2[20 * i]);
208     
209      diagptable = &diagptable_start[19 * cptr[i]];               
210     
211      term =  left[0] * right[0];           
212      term += left[1] * right[1] * *diagptable++; 
213      term += left[2] * right[2] * *diagptable++; 
214      term += left[3] * right[3] * *diagptable++; 
215      term += left[4] * right[4] * *diagptable++;         
216      term += left[5] * right[5] * *diagptable++; 
217      term += left[6] * right[6] * *diagptable++; 
218      term += left[7] * right[7] * *diagptable++; 
219      term += left[8] * right[8] * *diagptable++;
220      term += left[9] * right[9] * *diagptable++;           
221      term += left[10] * right[10] * *diagptable++; 
222      term += left[11] * right[11] * *diagptable++; 
223      term += left[12] * right[12] * *diagptable++; 
224      term += left[13] * right[13] * *diagptable++;
225      term += left[14] * right[14] * *diagptable++;         
226      term += left[15] * right[15] * *diagptable++; 
227      term += left[16] * right[16] * *diagptable++; 
228      term += left[17] * right[17] * *diagptable++; 
229      term += left[18] * right[18] * *diagptable++;
230      term += left[19] * right[19] * *diagptable++; 
231           
232      term = log(term) + ex2[i] * log(minlikelihood);             
233
234      v[i] = term;       
235    }
236
237  free(diagptable_start);         
238}
239     
240   
241static void evaluateGTRCATPROTMULT_VECTOR (int *ex2, int *cptr, int *modelptr, 
242                                           double *x2,  double *v, double *EIGN, double *rptr, double *tipVector, 
243                                           double *pz, 
244                                           char *tipX1, int lower, int n, int numberOfCategories, int numberOfModels, int multiBranch)
245{
246  double   z, lz = 0.0, ki, lza[19];       
247  int     i, l; 
248  double  *diagptable, *diagptable_start;   
249  double term;
250  double *left, *right;   
251  int model, modelCounter;     
252   
253  if(!multiBranch)
254    {
255      z = pz[0];   
256      if (z < zmin) z = zmin;
257      lz = log(z); 
258    }
259
260  diagptable = diagptable_start = (double *)malloc(sizeof(double) * numberOfCategories * 19 * numberOfModels); 
261 
262  for(modelCounter = 0; modelCounter < numberOfModels; modelCounter++)
263    {
264      if(!multiBranch)
265        {
266          z = pz[modelCounter];   
267          if (z < zmin) z = zmin;
268          lz = log(z); 
269        }
270     
271      for(l = 0; l < 19; l++)     
272         lza[l] = EIGN[modelCounter * 19 + l] * lz;       
273
274      for(i = 0; i <  numberOfCategories; i++)
275        {       
276          ki = rptr[i]; 
277         
278          *diagptable++ = exp (ki * lza[0]);       
279          *diagptable++ = exp (ki * lza[1]);
280          *diagptable++ = exp (ki * lza[2]);       
281          *diagptable++ = exp (ki * lza[3]);
282          *diagptable++ = exp (ki * lza[4]);       
283          *diagptable++ = exp (ki * lza[5]);
284          *diagptable++ = exp (ki * lza[6]);       
285          *diagptable++ = exp (ki * lza[7]);
286          *diagptable++ = exp (ki * lza[8]);       
287          *diagptable++ = exp (ki * lza[9]);
288          *diagptable++ = exp (ki * lza[10]);     
289          *diagptable++ = exp (ki * lza[11]);
290          *diagptable++ = exp (ki * lza[12]);     
291          *diagptable++ = exp (ki * lza[13]);
292          *diagptable++ = exp (ki * lza[14]);     
293          *diagptable++ = exp (ki * lza[15]);
294          *diagptable++ = exp (ki * lza[16]);     
295          *diagptable++ = exp (ki * lza[17]);
296          *diagptable++ = exp (ki * lza[18]);             
297        }
298    }                   
299               
300
301  for(i = lower; i < n; i++) 
302    {                       
303      model = modelptr[i];
304      left  = &(tipVector[model * 460 + 20 * tipX1[i]]);
305      right = &(x2[20 * i]);
306     
307      diagptable = &diagptable_start[model * 19 * numberOfCategories + 19 * cptr[i]];             
308     
309      term =  left[0] * right[0];           
310      term += left[1] * right[1] * *diagptable++; 
311      term += left[2] * right[2] * *diagptable++; 
312      term += left[3] * right[3] * *diagptable++; 
313      term += left[4] * right[4] * *diagptable++;         
314      term += left[5] * right[5] * *diagptable++; 
315      term += left[6] * right[6] * *diagptable++; 
316      term += left[7] * right[7] * *diagptable++; 
317      term += left[8] * right[8] * *diagptable++;
318      term += left[9] * right[9] * *diagptable++;           
319      term += left[10] * right[10] * *diagptable++; 
320      term += left[11] * right[11] * *diagptable++; 
321      term += left[12] * right[12] * *diagptable++; 
322      term += left[13] * right[13] * *diagptable++;
323      term += left[14] * right[14] * *diagptable++;         
324      term += left[15] * right[15] * *diagptable++; 
325      term += left[16] * right[16] * *diagptable++; 
326      term += left[17] * right[17] * *diagptable++; 
327      term += left[18] * right[18] * *diagptable++;
328      term += left[19] * right[19] * *diagptable++; 
329     
330      term = log(term) + ex2[i] * log(minlikelihood);             
331     
332      v[i] = term;         
333    }
334 
335  free(diagptable_start);             
336}         
337     
338static void evaluateGTRGAMMA_VECTOR (int *ex2,
339                                     double *x2_start, double *v, double *EIGN, double *gammaRates, 
340                                     double *tipVector, double pz, 
341                                     char *tipX1, int lower, int n)
342{
343  double   z, lz, term, ki;   
344  int     i; 
345  double  *diagptable, *diagptable_start;
346  double *x1, *x2; 
347 
348  z = pz;
349 
350  if (z < zmin) z = zmin;
351  lz = log(z);
352 
353  diagptable = diagptable_start = (double *)malloc(sizeof(double) * 4 * 3);
354
355  for(i = 0; i < 4; i++)
356    {
357      ki = gammaRates[i];       
358     
359      *diagptable++ = exp (EIGN[0] * ki * lz);
360      *diagptable++ = exp (EIGN[1] * ki * lz);
361      *diagptable++ = exp (EIGN[2] * ki * lz);
362    }                     
363 
364   
365  for (i = lower; i < n; i++) 
366    {
367      x1 = &(tipVector[4 * tipX1[i]]);
368      x2 = &x2_start[16 * i];
369      diagptable = diagptable_start;
370     
371      /* cat 0 */
372     
373      term =  x1[0] * x2[0];     
374      term += x1[1] * x2[1] * *diagptable++;   
375      term += x1[2] * x2[2] * *diagptable++;     
376      term += x1[3] * x2[3] * *diagptable++;     
377     
378      /* cat 1 */
379     
380      term += x1[0] * x2[4];
381      term += x1[1] * x2[5] * *diagptable++;
382      term += x1[2] * x2[6] * *diagptable++;
383      term += x1[3] * x2[7] * *diagptable++;     
384     
385      /* cat 2 */
386     
387      term += x1[0] * x2[8];
388      term += x1[1] * x2[9] * *diagptable++;
389      term += x1[2] * x2[10] * *diagptable++;
390      term += x1[3] * x2[11] * *diagptable++;     
391     
392      /* cat 3 */
393     
394      term += x1[0] * x2[12];
395      term += x1[1] * x2[13] * *diagptable++;
396      term += x1[2] * x2[14] * *diagptable++;
397      term += x1[3] * x2[15] * *diagptable++;     
398     
399      term = log(0.25 * term) + ex2[i] * log(minlikelihood);
400      v[i] = term;     
401    }
402 
403  free(diagptable_start);   
404} 
405
406static void evaluateGTRGAMMAMULT_VECTOR(int *ex2, int *modelptr, 
407                                        double *x2_start, double *v, double *EIGN, double *gammaRates, 
408                                        double *tipVector, double *pz, 
409                                        char *tipX1, int lower, int n, int numberOfModels, int multiBranch)
410{
411  double   z, lz = 0.0, term, ki;   
412  int     i;
413  double  *diagptable, *diagptable_start; 
414  double *x1, *x2;
415  int model;
416 
417  if(!multiBranch)
418    {
419      z = pz[0]; 
420      if (z < zmin) z = zmin;
421      lz = log(z);
422    }
423 
424  diagptable = diagptable_start = (double *)malloc(sizeof(double) * 12 * numberOfModels); 
425
426  for(model = 0; model < numberOfModels; model++)
427    {
428      if(multiBranch)
429        {
430          z = pz[model]; 
431          if (z < zmin) z = zmin;
432          lz = log(z);
433        }
434     
435      diagptable = &diagptable_start[12 * model];
436     
437      for(i = 0; i < 4; i++)
438        {
439          ki = gammaRates[model * 4 + i];       
440         
441          *diagptable++ = exp (EIGN[model * 3] * ki * lz);
442          *diagptable++ = exp (EIGN[model * 3 + 1] * ki * lz);
443          *diagptable++ = exp (EIGN[model * 3 + 2] * ki * lz);
444        }
445    }               
446
447
448  for (i = lower; i < n; i++) 
449    {
450      model = modelptr[i];
451     
452      x1 = &(tipVector[64 * model + 4 * tipX1[i]]);     
453      x2 = &x2_start[16 * i];
454     
455      diagptable = &diagptable_start[12 * model];
456     
457      /* cat 0 */
458     
459      term =  x1[0] * x2[0];     
460      term += x1[1] * x2[1] * *diagptable++;   
461      term += x1[2] * x2[2] * *diagptable++;     
462      term += x1[3] * x2[3] * *diagptable++;     
463     
464      /* cat 1 */
465     
466      term += x1[0] * x2[4];
467      term += x1[1] * x2[5] * *diagptable++;
468      term += x1[2] * x2[6] * *diagptable++;
469      term += x1[3] * x2[7] * *diagptable++;     
470     
471      /* cat 2 */
472     
473      term += x1[0] * x2[8];
474      term += x1[1] * x2[9] * *diagptable++;
475      term += x1[2] * x2[10] * *diagptable++;
476      term += x1[3] * x2[11] * *diagptable++;     
477     
478      /* cat 3 */
479     
480      term += x1[0] * x2[12];
481      term += x1[1] * x2[13] * *diagptable++;
482      term += x1[2] * x2[14] * *diagptable++;
483      term += x1[3] * x2[15] * *diagptable++;     
484     
485      term = log(0.25 * term) + ex2[i] * log(minlikelihood);
486      v[i] = term;         
487    }
488         
489  free(diagptable_start); 
490}
491     
492static void evaluateGTRGAMMAMULTINVAR_VECTOR(int *ex2, int *modelptr, int *iptr,
493                                             double *x2_start, double *v, double *EIGN, double *gammaRates, 
494                                             double *tipVector, double *tFreqs, double *invariants, double *pz, 
495                                             char *tipX1, int lower, int n, int numberOfModels, int multiBranch)
496{
497  double   z, lz = 0.0, term, ki;   
498  int     i, model; 
499  double  *diagptable, *diagptable_start;
500  double *x1, *x2; 
501  double *scalers;
502  double *freqs;   
503 
504  if(!multiBranch)
505    {
506      z = pz[0]; 
507      if (z < zmin) z = zmin;
508      lz = log(z);
509    }
510 
511  diagptable = diagptable_start = (double *)malloc(sizeof(double) * 12 * numberOfModels);
512  freqs      = (double *)malloc(4 * numberOfModels * sizeof(double));
513  scalers    = (double *)malloc(numberOfModels * sizeof(double)); 
514
515  for(model = 0; model < numberOfModels; model++)
516    {
517      if(multiBranch)
518        {
519          z = pz[model]; 
520          if (z < zmin) z = zmin;
521          lz = log(z);
522        }
523 
524      scalers[model] = 0.25 * (1.0 - invariants[model]); 
525      freqs[4 * model]     = tFreqs[4 * model]     * invariants[model]; 
526      freqs[4 * model + 1] = tFreqs[4 * model + 1] * invariants[model];
527      freqs[4 * model + 2] = tFreqs[4 * model + 2] * invariants[model];
528      freqs[4 * model + 3] = tFreqs[4 * model + 3] * invariants[model];
529     
530      diagptable = &diagptable_start[12 * model];
531     
532      for(i = 0; i < 4; i++)
533        {
534          ki = gammaRates[model * 4 + i];       
535         
536          *diagptable++ = exp (EIGN[model * 3] * ki * lz);
537          *diagptable++ = exp (EIGN[model * 3 + 1] * ki * lz);
538          *diagptable++ = exp (EIGN[model * 3 + 2] * ki * lz);
539        }
540    }
541               
542
543  for (i = lower; i < n; i++) 
544    {
545      model = modelptr[i];
546     
547      x1 = &(tipVector[64 * model + 4 * tipX1[i]]);     
548      x2 = &x2_start[16 * i];
549     
550      diagptable = &diagptable_start[12 * model];
551     
552      /* cat 0 */
553     
554      term =  x1[0] * x2[0];     
555      term += x1[1] * x2[1] * *diagptable++;   
556      term += x1[2] * x2[2] * *diagptable++;     
557      term += x1[3] * x2[3] * *diagptable++;     
558     
559      /* cat 1 */
560     
561      term += x1[0] * x2[4];
562      term += x1[1] * x2[5] * *diagptable++;
563      term += x1[2] * x2[6] * *diagptable++;
564      term += x1[3] * x2[7] * *diagptable++;     
565     
566      /* cat 2 */
567     
568      term += x1[0] * x2[8];
569      term += x1[1] * x2[9] * *diagptable++;
570      term += x1[2] * x2[10] * *diagptable++;
571      term += x1[3] * x2[11] * *diagptable++;     
572     
573      /* cat 3 */
574     
575      term += x1[0] * x2[12];
576      term += x1[1] * x2[13] * *diagptable++;
577      term += x1[2] * x2[14] * *diagptable++;
578      term += x1[3] * x2[15] * *diagptable++;             
579     
580      if(iptr[i] < 4)     
581        term = log((scalers[model] * term + freqs[model * 4 + iptr[i]]) * pow(minlikelihood, ex2[i]));
582      else
583        term = log(scalers[model] * term) + ex2[i] * log(minlikelihood);
584       
585      v[i] = term;     
586    }     
587         
588  free(diagptable_start); 
589  free(scalers);
590  free(freqs);
591}
592     
593static void evaluateGTRGAMMAINVAR_VECTOR(int *ex2, int *iptr,
594                                         double *x2_start, double *v, double *EIGN, double *gammaRates, 
595                                         double *tipVector, double *tFreqs, double invariants, double pz, 
596                                         char *tipX1, int lower, int n)
597{
598  double   z, lz, term, ki;   
599  int     i; 
600  double  *diagptable; 
601  double freqs[4]; 
602  double *x1, *x2;
603  double scaler = 0.25 * (1.0 - invariants);   
604
605  freqs[0] = tFreqs[0] * invariants; 
606  freqs[1] = tFreqs[1] * invariants;
607  freqs[2] = tFreqs[2] * invariants;
608  freqs[3] = tFreqs[3] * invariants;
609 
610  z = pz;
611 
612  if (z < zmin) z = zmin;
613 
614  lz = log(z);
615 
616  diagptable = (double *)malloc(sizeof(double) * 16); 
617
618  for(i = 0; i < 4; i++)
619    {
620      ki = gammaRates[i];       
621     
622      diagptable[i * 4]     = 1.0;
623      diagptable[i * 4 + 1] = exp (EIGN[0] * ki * lz);
624      diagptable[i * 4 + 2] = exp (EIGN[1] * ki * lz);
625      diagptable[i * 4 + 3] = exp (EIGN[2] * ki * lz);
626    }                 
627
628
629  for (i = lower; i < n; i++) 
630    {
631      x1 = &(tipVector[4 * tipX1[i]]);
632      x2 = &x2_start[16 * i];
633     
634      /* cat 0 */
635     
636      term =  x1[0] * x2[0] * diagptable[0];     
637      term += x1[1] * x2[1] * diagptable[1];   
638      term += x1[2] * x2[2] * diagptable[2];     
639      term += x1[3] * x2[3] * diagptable[3];     
640     
641      /* cat 1 */
642     
643      term += x1[0] * x2[4] * diagptable[4];
644      term += x1[1] * x2[5] * diagptable[5];
645      term += x1[2] * x2[6] * diagptable[6];
646      term += x1[3] * x2[7] * diagptable[7];     
647     
648      /* cat 2 */
649     
650      term += x1[0] * x2[8] * diagptable[8];
651      term += x1[1] * x2[9] * diagptable[9];
652      term += x1[2] * x2[10] * diagptable[10];
653      term += x1[3] * x2[11] * diagptable[11];     
654     
655      /* cat 3 */
656     
657      term += x1[0] * x2[12] * diagptable[12];
658      term += x1[1] * x2[13] * diagptable[13];
659      term += x1[2] * x2[14] * diagptable[14];
660      term += x1[3] * x2[15] * diagptable[15];   
661     
662      if(iptr[i] < 4)
663        term = log((scaler * term + freqs[iptr[i]]) * pow(minlikelihood, ex2[i]));
664      else
665        term = log(scaler * term) + ex2[i] * log(minlikelihood);         
666     
667      v[i] = term;                             
668    }
669         
670  free(diagptable);           
671}
672 
673static void evaluateGTRGAMMAPROT_VECTOR(int *ex2, 
674                                        double *x2,  double *v, double *EIGN, double *gammaRates, 
675                                        double *tipVector, double pz, 
676                                        char *tipX1, int lower, int n)
677{
678  double   z, lz, term, ki;   
679  int     i, j;
680  double  *diagptable, *diagptable_start;
681  double  *left, *right; 
682         
683  z = pz;
684   
685  if (z < zmin) z = zmin;
686  lz = log(z);
687 
688  diagptable = diagptable_start = (double *)malloc(sizeof(double) * 4 * 19);
689
690  for(i = 0; i < 4; i++)
691    {
692      ki = gammaRates[i];       
693     
694      for(j = 0; j < 19; j++)
695        *diagptable++ = exp (EIGN[j] * ki * lz);       
696    }                       
697
698
699  for (i = lower; i < n; i++) 
700    {
701      left = &(tipVector[20 * tipX1[i]]);
702      diagptable = diagptable_start;
703     
704      term = 0;
705     
706      for(j = 0; j < 4; j++)
707        {
708          right = &(x2[80 * i + 20 * j]);
709         
710          term +=  left[0] * right[0];
711          term +=  left[1] * right[1] * *diagptable++;
712          term +=  left[2] * right[2] * *diagptable++;
713          term +=  left[3] * right[3] * *diagptable++;
714          term +=  left[4] * right[4] * *diagptable++;
715          term +=  left[5] * right[5] * *diagptable++;
716          term +=  left[6] * right[6] * *diagptable++; 
717          term +=  left[7] * right[7] * *diagptable++;
718          term +=  left[8] * right[8] * *diagptable++;
719          term +=  left[9] * right[9] * *diagptable++;
720          term +=  left[10] * right[10] * *diagptable++;
721          term +=  left[11] * right[11] * *diagptable++;
722          term +=  left[12] * right[12] * *diagptable++;
723          term +=  left[13] * right[13] * *diagptable++;
724          term +=  left[14] * right[14] * *diagptable++;
725          term +=  left[15] * right[15] * *diagptable++;
726          term +=  left[16] * right[16] * *diagptable++;
727          term +=  left[17] * right[17] * *diagptable++;
728          term +=  left[18] * right[18] * *diagptable++;       
729          term +=  left[19] * right[19] * *diagptable++;
730        }         
731     
732      term = log(0.25 * term) + ex2[i] * log(minlikelihood);       
733      v[i] = term;                               
734    }
735   
736  free(diagptable_start); 
737}
738   
739static void evaluateGTRGAMMAPROTINVAR_VECTOR(int *ex2, int *iptr,
740                                             double *x2,  double *v, double *EIGN, double *gammaRates, 
741                                             double *tipVector,double *tFreqs, double invariants, double pz, 
742                                             char *tipX1, int lower, int n)
743{
744  double   z, lz, term, ki;   
745  int     i, j;   
746  double  *diagptable, *diagptable_start;
747  double *left, *right;
748  double scaler = 0.25 * (1.0 - invariants);
749  double freqs[20]; 
750 
751  for(i = 0; i < 20; i++)
752    freqs[i] = tFreqs[i] * invariants;
753   
754  z = pz;
755 
756  if (z < zmin) z = zmin;
757  lz = log(z);
758 
759  diagptable = diagptable_start = (double *)malloc(sizeof(double) * 4 * 19);
760
761  for(i = 0; i < 4; i++)
762    {
763      ki = gammaRates[i];       
764     
765      for(j = 0; j < 19; j++)
766        *diagptable++ = exp (EIGN[j] * ki * lz);       
767    }               
768
769
770  for (i = lower; i < n; i++) 
771    {
772      left = &(tipVector[20 * tipX1[i]]);
773      diagptable = diagptable_start;
774     
775      term = 0;
776
777      for(j = 0; j < 4; j++)
778        {
779          right = &(x2[80 * i + 20 * j]);
780         
781          term +=  left[0] * right[0];
782          term +=  left[1] * right[1] * *diagptable++;
783          term +=  left[2] * right[2] * *diagptable++;
784          term +=  left[3] * right[3] * *diagptable++;
785          term +=  left[4] * right[4] * *diagptable++;
786          term +=  left[5] * right[5] * *diagptable++;
787          term +=  left[6] * right[6] * *diagptable++; 
788          term +=  left[7] * right[7] * *diagptable++;
789          term +=  left[8] * right[8] * *diagptable++;
790          term +=  left[9] * right[9] * *diagptable++;
791          term +=  left[10] * right[10] * *diagptable++;
792          term +=  left[11] * right[11] * *diagptable++;
793          term +=  left[12] * right[12] * *diagptable++;
794          term +=  left[13] * right[13] * *diagptable++;
795          term +=  left[14] * right[14] * *diagptable++;
796          term +=  left[15] * right[15] * *diagptable++;
797          term +=  left[16] * right[16] * *diagptable++;
798          term +=  left[17] * right[17] * *diagptable++;
799          term +=  left[18] * right[18] * *diagptable++;       
800          term +=  left[19] * right[19] * *diagptable++;
801        }         
802     
803      if(iptr[i] < 20)     
804        term = log((scaler * term + freqs[iptr[i]]) * pow(minlikelihood, ex2[i]));
805      else
806        term = log(scaler * term) + ex2[i] * log(minlikelihood);
807     
808      v[i] = term;                               
809    }
810   
811
812  free(diagptable_start); 
813}
814         
815static void evaluateGTRGAMMAPROTMULT_VECTOR (int *ex2, int *modelptr,
816                                             double *x2,  double *v, double *EIGN, double *gammaRates, 
817                                             double *tipVector, double *pz, 
818                                             char *tipX1, int lower, int n, int numberOfModels, int multiBranch)
819{
820  double   z, lz = 0.0, term, ki;   
821  int     i, j; 
822  double  *diagptable, *diagptable_start;
823  double *left, *right;
824  int model;
825     
826  if(!multiBranch)
827    {
828      z = pz[0]; 
829      if (z < zmin) z = zmin;
830      lz = log(z);
831    }
832 
833  diagptable = diagptable_start = (double *)malloc(sizeof(double) * 4 * 19 * numberOfModels);
834 
835  for(model = 0; model < numberOfModels; model++)
836    {
837       if(multiBranch)
838         {
839           z = pz[model]; 
840           if (z < zmin) z = zmin;
841           lz = log(z);
842         }
843     
844      diagptable = &diagptable_start[4 * 19 * model];
845     
846      for(i = 0; i < 4; i++)
847        {
848          ki = gammaRates[model * 4 + i];       
849         
850          for(j = 0; j < 19; j++)
851            *diagptable++ = exp (EIGN[model * 19 + j] * ki * lz);       
852        }
853    }                           
854             
855
856  for (i = lower; i < n; i++) 
857    {       
858      model = modelptr[i];
859
860      left = &(tipVector[460 * model + 20 * tipX1[i]]);
861      diagptable = &diagptable_start[76 * model];         
862     
863      term = 0;
864     
865      for(j = 0; j < 4; j++)
866        {
867          right = &(x2[80 * i + 20 * j]);
868         
869          term +=  left[0] * right[0];
870          term +=  left[1] * right[1] * *diagptable++;
871          term +=  left[2] * right[2] * *diagptable++;
872          term +=  left[3] * right[3] * *diagptable++;
873          term +=  left[4] * right[4] * *diagptable++;
874          term +=  left[5] * right[5] * *diagptable++;
875          term +=  left[6] * right[6] * *diagptable++; 
876          term +=  left[7] * right[7] * *diagptable++;
877          term +=  left[8] * right[8] * *diagptable++;
878          term +=  left[9] * right[9] * *diagptable++;
879          term +=  left[10] * right[10] * *diagptable++;
880          term +=  left[11] * right[11] * *diagptable++;
881          term +=  left[12] * right[12] * *diagptable++;
882          term +=  left[13] * right[13] * *diagptable++;
883          term +=  left[14] * right[14] * *diagptable++;
884          term +=  left[15] * right[15] * *diagptable++;
885          term +=  left[16] * right[16] * *diagptable++;
886          term +=  left[17] * right[17] * *diagptable++;
887          term +=  left[18] * right[18] * *diagptable++;       
888          term +=  left[19] * right[19] * *diagptable++;
889        }         
890     
891      term = log(0.25 * term) + ex2[i] * log(minlikelihood);       
892      v[i] = term;                               
893    }
894   
895  free(diagptable_start); 
896}
897       
898static void evaluateGTRGAMMAPROTMULTINVAR_VECTOR(int *ex2, int *modelptr, int *iptr,
899                                                 double *x2,  double *v, double *EIGN, double *gammaRates, 
900                                                 double *tipVector, double *tFreqs, double *invariants, double *pz, 
901                                                 char *tipX1, int lower, int n, int numberOfModels, int multiBranch)
902{
903  double   z, lz = 0.0, term, ki;   
904  int     i, j; 
905  double  *diagptable, *diagptable_start;
906  double *left, *right;
907  int model;
908  double *scalers;
909  double *freqs; 
910     
911  if(!multiBranch)
912    {
913      z = pz[0]; 
914      if (z < zmin) z = zmin;
915      lz = log(z);
916    }
917 
918  diagptable = diagptable_start = (double *)malloc(sizeof(double) * 4 * 19 * numberOfModels);
919  freqs      = (double *)malloc(20 * numberOfModels * sizeof(double));
920  scalers    = (double *)malloc(numberOfModels * sizeof(double));
921 
922  for(model = 0; model < numberOfModels; model++)
923    {
924      if(multiBranch)
925        {
926          z = pz[model]; 
927          if (z < zmin) z = zmin;
928          lz = log(z);
929        }       
930           
931      diagptable = &diagptable_start[4 * 19 * model];
932
933      scalers[model] = 0.25 * (1.0 - invariants[model]); 
934     
935      for(i = 0; i < 20; i++)
936        freqs[20 * model + i] = tFreqs[20 * model + i] * invariants[model]; 
937     
938      for(i = 0; i < 4; i++)
939        {
940          ki = gammaRates[model * 4 + i];       
941         
942          for(j = 0; j < 19; j++)
943            *diagptable++ = exp (EIGN[model * 19 + j] * ki * lz);       
944        }
945    }                 
946
947       
948  for (i = lower; i < n; i++) 
949    {       
950      model = modelptr[i];
951     
952      left = &(tipVector[460 * model + 20 * tipX1[i]]);
953      diagptable = &diagptable_start[76 * model];         
954     
955      term = 0;
956     
957      for(j = 0; j < 4; j++)
958        {
959          right = &(x2[80 * i + 20 * j]);
960         
961          term +=  left[0] * right[0];
962          term +=  left[1] * right[1] * *diagptable++;
963          term +=  left[2] * right[2] * *diagptable++;
964          term +=  left[3] * right[3] * *diagptable++;
965          term +=  left[4] * right[4] * *diagptable++;
966          term +=  left[5] * right[5] * *diagptable++;
967          term +=  left[6] * right[6] * *diagptable++; 
968          term +=  left[7] * right[7] * *diagptable++;
969          term +=  left[8] * right[8] * *diagptable++;
970          term +=  left[9] * right[9] * *diagptable++;
971          term +=  left[10] * right[10] * *diagptable++;
972          term +=  left[11] * right[11] * *diagptable++;
973          term +=  left[12] * right[12] * *diagptable++;
974          term +=  left[13] * right[13] * *diagptable++;
975          term +=  left[14] * right[14] * *diagptable++;
976          term +=  left[15] * right[15] * *diagptable++;
977          term +=  left[16] * right[16] * *diagptable++;
978          term +=  left[17] * right[17] * *diagptable++;
979          term +=  left[18] * right[18] * *diagptable++;       
980          term +=  left[19] * right[19] * *diagptable++;
981        }                                         
982     
983      if(iptr[i] < 20)   
984        term = log((scalers[model] * term + freqs[model * 20 + iptr[i]]) * pow(minlikelihood, ex2[i]));
985      else
986        term = log(scalers[model] * term) + ex2[i] * log(minlikelihood);
987     
988      v[i] = term;                               
989    }
990   
991
992  free(diagptable_start); 
993  free(freqs);
994  free(scalers);
995}
996       
997   
998
999#ifdef _LOCAL_DATA
1000 
1001void evaluateGenericVectorIterative(tree *localTree, int startIndex, int endIndex)   
1002{
1003  double *x2_start, *pz, *v;
1004  char *tip;
1005  int *ex2;
1006  int pNumber, qNumber;
1007
1008  v = localTree->strided_siteLL_Vector;
1009
1010  pNumber = localTree->td[0].ti[0].pNumber;
1011  qNumber = localTree->td[0].ti[0].qNumber;
1012  pz      = localTree->td[0].ti[0].qz;
1013
1014  newviewIterative(localTree, startIndex, endIndex);
1015
1016
1017  x2_start = getLikelihoodArray(qNumber,  localTree->mxtips, localTree->xVector); 
1018  ex2      = getScalingArray(qNumber, localTree->cdta->endsite, localTree->mxtips, localTree->expArray); 
1019
1020  tip   = localTree->strided_yVector[pNumber];
1021
1022
1023  if(localTree->mixedData)
1024    {
1025      assert(0);
1026    }
1027  else
1028    {
1029      switch(localTree->likelihoodFunction)
1030        {
1031        case GTRCAT:
1032          evaluateGTRCAT_VECTOR(ex2, localTree->strided_rateCategory, 
1033                                x2_start, v, localTree->EIGN_DNA, localTree->strided_patrat, localTree->tipVectorDNA, pz[0], 
1034                                tip, startIndex, endIndex, localTree->NumberOfCategories);   
1035          break;
1036        case GTRCATMULT:
1037          evaluateGTRCATMULT_VECTOR(ex2, localTree->strided_rateCategory, localTree->strided_model,
1038                                    x2_start, v, localTree->EIGN_DNA, localTree->strided_patrat, localTree->tipVectorDNA, pz, 
1039                                    tip, startIndex, endIndex, localTree->NumberOfCategories, localTree->NumberOfModels, 
1040                                    localTree->multiBranch);
1041          break;
1042        case PROTCAT:
1043          evaluateGTRCATPROT_VECTOR(ex2, localTree->strided_rateCategory,
1044                                    x2_start, v, localTree->EIGN_AA, localTree->strided_patrat, localTree->tipVectorAA, pz[0], 
1045                                    tip, startIndex, endIndex, localTree->NumberOfCategories);
1046          break;
1047        case PROTCATMULT:
1048          evaluateGTRCATPROTMULT_VECTOR(ex2, localTree->strided_rateCategory, localTree->strided_model, 
1049                                        x2_start, v, localTree->EIGN_AA, localTree->strided_patrat, localTree->tipVectorAA, pz, 
1050                                        tip, startIndex, endIndex, localTree->NumberOfCategories, localTree->NumberOfModels, localTree->multiBranch);
1051          break;
1052        case GTRGAMMA:
1053          evaluateGTRGAMMA_VECTOR(ex2, 
1054                                  x2_start, v, localTree->EIGN_DNA, localTree->gammaRates, localTree->tipVectorDNA, pz[0], 
1055                                  tip, startIndex, endIndex); 
1056          break;
1057        case GTRGAMMAI:
1058          evaluateGTRGAMMAINVAR_VECTOR(ex2, localTree->strided_invariant,
1059                                       x2_start, v, localTree->EIGN_DNA, localTree->gammaRates, 
1060                                       localTree->tipVectorDNA, localTree->frequencies_DNA, localTree->invariants[0], pz[0], 
1061                                       tip, startIndex, endIndex);
1062          break;
1063        case GTRGAMMAMULT:
1064          evaluateGTRGAMMAMULT_VECTOR(ex2, localTree->strided_model,
1065                                      x2_start, v, localTree->EIGN_DNA, localTree->gammaRates, localTree->tipVectorDNA, pz, 
1066                                      tip, startIndex, endIndex, localTree->NumberOfModels, localTree->multiBranch);
1067          break;
1068        case GTRGAMMAMULTI:
1069          evaluateGTRGAMMAMULTINVAR_VECTOR(ex2, localTree->strided_model,  localTree->strided_invariant,
1070                                           x2_start, v, localTree->EIGN_DNA, localTree->gammaRates, 
1071                                           localTree->tipVectorDNA, localTree->frequencies_DNA, localTree->invariants, pz, 
1072                                           tip, startIndex, endIndex, localTree->NumberOfModels, localTree->multiBranch); 
1073          break;
1074        case PROTGAMMA:
1075          evaluateGTRGAMMAPROT_VECTOR(ex2,
1076                                      x2_start, v, localTree->EIGN_AA, localTree->gammaRates, localTree->tipVectorAA, pz[0], 
1077                                      tip, startIndex, endIndex);
1078          break;
1079        case PROTGAMMAI:
1080          evaluateGTRGAMMAPROTINVAR_VECTOR(ex2,  localTree->strided_invariant,
1081                                           x2_start, v, localTree->EIGN_AA, localTree->gammaRates, 
1082                                           localTree->tipVectorAA, localTree->frequencies_AA, localTree->invariants[0], pz[0], 
1083                                           tip, startIndex, endIndex);
1084          break;
1085        case PROTGAMMAMULT:
1086          evaluateGTRGAMMAPROTMULT_VECTOR(ex2, localTree->strided_model,
1087                                          x2_start, v, localTree->EIGN_AA, localTree->gammaRates, localTree->tipVectorAA, pz, 
1088                                          tip, startIndex, endIndex, localTree->NumberOfModels, localTree->multiBranch);
1089          break;
1090        case PROTGAMMAMULTI:
1091          evaluateGTRGAMMAPROTMULTINVAR_VECTOR(ex2, localTree->strided_model,  localTree->strided_invariant,
1092                                               x2_start, v, localTree->EIGN_AA, localTree->gammaRates, 
1093                                               localTree->tipVectorAA, localTree->frequencies_AA, localTree->invariants, pz, 
1094                                               tip, startIndex, endIndex, localTree->NumberOfModels, localTree->multiBranch
1095                                               );
1096          break;
1097        default:
1098          assert(0);
1099        } 
1100    }
1101
1102}
1103 
1104#else
1105
1106void evaluateGenericVectorIterative(tree *tr, int startIndex, int endIndex)   
1107{
1108  double *x2_start, *pz, *v;
1109  char *tip;
1110  int *ex2;
1111  int pNumber, qNumber;
1112
1113  v = tr->siteLL_Vector;
1114
1115  pNumber = tr->td[0].ti[0].pNumber;
1116  qNumber = tr->td[0].ti[0].qNumber;
1117  pz      = tr->td[0].ti[0].qz;
1118
1119  newviewIterative(tr, 0, tr->cdta->endsite);
1120
1121
1122  x2_start = getLikelihoodArray(qNumber,  tr->mxtips, tr->xVector); 
1123  ex2      = getScalingArray(qNumber, tr->cdta->endsite, tr->mxtips, tr->expArray); 
1124
1125  tip   = tr->yVector[pNumber];
1126
1127
1128  if(tr->mixedData)
1129    {
1130      int model, branchIndex, width, l, u, offset;
1131
1132      for(model = 0; model < tr->NumberOfModels; model++)
1133        {
1134          if(tr->multiBranch)
1135            branchIndex = model;
1136          else
1137            branchIndex = 0;
1138
1139          l = tr->partitionData[model].lower;
1140          u = tr->partitionData[model].upper;
1141         
1142          offset =  tr->partitionData[model].modelOffset;
1143
1144          width = u - l;
1145
1146          switch(tr->partitionData[model].dataType)
1147            {
1148            case DNA_DATA:
1149              switch(tr->rateHetModel)
1150                {
1151                case CAT:
1152                  evaluateGTRCAT_VECTOR(&(ex2[l]), &(tr->cdta->rateCategory[l]), 
1153                                        &(x2_start[offset]), &(v[l]), &(tr->EIGN_DNA[model * 3]), 
1154                                        tr->cdta->patrat, &(tr->tipVectorDNA[64 * model]), pz[branchIndex], 
1155                                        &(tip[l]), 0, width, tr->NumberOfCategories);
1156                  break;
1157                case GAMMA:               
1158                  evaluateGTRGAMMA_VECTOR(&(ex2[l]), 
1159                                          &(x2_start[offset]), &(v[l]), 
1160                                          &(tr->EIGN_DNA[model * 3]), 
1161                                          &(tr->gammaRates[model * 4]), 
1162                                          &(tr->tipVectorDNA[model * 64]), pz[branchIndex], 
1163                                          &(tip[l]), 0, width); 
1164                  break;
1165                case GAMMA_I:
1166                  evaluateGTRGAMMAINVAR_VECTOR(&(ex2[l]), &(tr->invariant[l]),
1167                                               &(x2_start[offset]), &(v[l]), 
1168                                               &(tr->EIGN_DNA[model * 3]), 
1169                                               &(tr->gammaRates[model * 4]), 
1170                                               &(tr->tipVectorDNA[model * 64]), 
1171                                               &(tr->frequencies_DNA[model * 4]), 
1172                                               tr->invariants[model], pz[branchIndex], 
1173                                               &(tip[l]), 0, width);
1174                  break;
1175                default:
1176                  assert(0);
1177                }
1178              break;
1179            case AA_DATA:
1180              switch(tr->rateHetModel)
1181                {
1182                case CAT:
1183                  evaluateGTRCATPROT_VECTOR(&(ex2[l]), &(tr->cdta->rateCategory[l]),
1184                                            &(x2_start[offset]), &(v[l]),
1185                                            &(tr->EIGN_AA[model * 19]), 
1186                                            tr->cdta->patrat, &(tr->tipVectorAA[model * 460]), pz[branchIndex], 
1187                                            &(tip[l]), 0, width, tr->NumberOfCategories);
1188                  break;
1189                case GAMMA:               
1190                  evaluateGTRGAMMAPROT_VECTOR(&(ex2[l]),
1191                                              &(x2_start[offset]), &(v[l]), 
1192                                              &(tr->EIGN_AA[model * 19]), 
1193                                              &(tr->gammaRates[model * 4]), 
1194                                              &(tr->tipVectorAA[model * 460]), pz[branchIndex], 
1195                                              &(tip[l]), 0, width);
1196                  break;
1197                case GAMMA_I:           
1198                  evaluateGTRGAMMAPROTINVAR_VECTOR(&(ex2[l]),  &(tr->invariant[l]),
1199                                                   &(x2_start[offset]), &(v[l]), 
1200                                                   &(tr->EIGN_AA[model * 19]), 
1201                                                   &(tr->gammaRates[model * 4]), 
1202                                                   &(tr->tipVectorAA[model * 460]), 
1203                                                   &(tr->frequencies_AA[model * 20]) , 
1204                                                   tr->invariants[model], pz[branchIndex], 
1205                                                   &(tip[l]), 0, width);
1206                  break;
1207                default:
1208                  assert(0);
1209                }
1210              break;
1211            default:
1212              assert(0);
1213            }
1214        }
1215     
1216    }
1217  else
1218    {
1219      switch(tr->likelihoodFunction)
1220        {
1221        case GTRCAT:
1222          evaluateGTRCAT_VECTOR(ex2, tr->cdta->rateCategory, 
1223                                x2_start, v, tr->EIGN_DNA, tr->cdta->patrat, tr->tipVectorDNA, pz[0], 
1224                                tip, startIndex, endIndex, tr->NumberOfCategories);   
1225          break;
1226        case GTRCATMULT:
1227          evaluateGTRCATMULT_VECTOR(ex2, tr->cdta->rateCategory, tr->model,
1228                                    x2_start, v, tr->EIGN_DNA, tr->cdta->patrat, tr->tipVectorDNA, pz, 
1229                                    tip, startIndex, endIndex, tr->NumberOfCategories, tr->NumberOfModels, tr->multiBranch);
1230          break;
1231        case PROTCAT:
1232          evaluateGTRCATPROT_VECTOR(ex2, tr->cdta->rateCategory,
1233                                    x2_start, v, tr->EIGN_AA, tr->cdta->patrat, tr->tipVectorAA, pz[0], 
1234                                    tip, startIndex, endIndex, tr->NumberOfCategories);
1235          break;
1236        case PROTCATMULT:
1237          evaluateGTRCATPROTMULT_VECTOR(ex2, tr->cdta->rateCategory, tr->model, 
1238                                        x2_start, v, tr->EIGN_AA, tr->cdta->patrat, tr->tipVectorAA, pz, 
1239                                        tip, startIndex, endIndex, tr->NumberOfCategories, tr->NumberOfModels, tr->multiBranch);
1240          break;
1241        case GTRGAMMA:
1242          evaluateGTRGAMMA_VECTOR(ex2, 
1243                                  x2_start, v, tr->EIGN_DNA, tr->gammaRates, tr->tipVectorDNA, pz[0], 
1244                                  tip, startIndex, endIndex); 
1245          break;
1246        case GTRGAMMAI:
1247          evaluateGTRGAMMAINVAR_VECTOR(ex2, tr->invariant,
1248                                       x2_start, v, tr->EIGN_DNA, tr->gammaRates, 
1249                                       tr->tipVectorDNA, tr->frequencies_DNA, tr->invariants[0], pz[0], 
1250                                       tip, startIndex, endIndex);
1251          break;
1252        case GTRGAMMAMULT:
1253          evaluateGTRGAMMAMULT_VECTOR(ex2, tr->model,
1254                                      x2_start, v, tr->EIGN_DNA, tr->gammaRates, tr->tipVectorDNA, pz, 
1255                                      tip, startIndex, endIndex, tr->NumberOfModels, tr->multiBranch);
1256          break;
1257        case GTRGAMMAMULTI:
1258          evaluateGTRGAMMAMULTINVAR_VECTOR(ex2, tr->model,  tr->invariant,
1259                                           x2_start, v, tr->EIGN_DNA, tr->gammaRates, 
1260                                           tr->tipVectorDNA, tr->frequencies_DNA, tr->invariants, pz, 
1261                                           tip, startIndex, endIndex, tr->NumberOfModels, tr->multiBranch); 
1262          break;
1263        case PROTGAMMA:
1264          evaluateGTRGAMMAPROT_VECTOR(ex2,
1265                                      x2_start, v, tr->EIGN_AA, tr->gammaRates, tr->tipVectorAA, pz[0], 
1266                                      tip, startIndex, endIndex);
1267          break;
1268        case PROTGAMMAI:
1269          evaluateGTRGAMMAPROTINVAR_VECTOR(ex2,  tr->invariant,
1270                                           x2_start, v, tr->EIGN_AA, tr->gammaRates, 
1271                                           tr->tipVectorAA, tr->frequencies_AA, tr->invariants[0], pz[0], 
1272                                           tip, startIndex, endIndex);
1273          break;
1274        case PROTGAMMAMULT:
1275          evaluateGTRGAMMAPROTMULT_VECTOR(ex2, tr->model,
1276                                          x2_start, v, tr->EIGN_AA, tr->gammaRates, tr->tipVectorAA, pz, 
1277                                          tip, startIndex, endIndex, tr->NumberOfModels, tr->multiBranch);
1278          break;
1279        case PROTGAMMAMULTI:
1280          evaluateGTRGAMMAPROTMULTINVAR_VECTOR(ex2, tr->model,  tr->invariant,
1281                                               x2_start, v, tr->EIGN_AA, tr->gammaRates, 
1282                                               tr->tipVectorAA, tr->frequencies_AA, tr->invariants, pz, 
1283                                               tip, startIndex, endIndex, tr->NumberOfModels, tr->multiBranch
1284                                               );
1285          break;
1286        default:
1287          assert(0);
1288        } 
1289    }
1290
1291}
1292
1293#endif
1294
1295void evaluateGenericVector (tree *tr, nodeptr p, double *v)
1296{
1297  nodeptr q = p->back;
1298  int i; 
1299   
1300  assert(isTip(p->number, tr->rdta->numsp)  || isTip(q->number, tr->rdta->numsp));
1301     
1302  if(isTip(q->number, tr->rdta->numsp))
1303    {
1304      nodeptr tmp = p;
1305      p = q;
1306      q = tmp;       
1307    }   
1308 
1309  tr->td[0].ti[0].pNumber = p->number;
1310  tr->td[0].ti[0].qNumber = q->number;
1311
1312  for(i = 0; i < tr->numBranches; i++)   
1313    tr->td[0].ti[0].qz[i] =  q->z[i];
1314
1315  tr->siteLL_Vector = v;
1316
1317  tr->td[0].count = 1;
1318  if(!q->x)
1319    computeTraversalInfo(q, &(tr->td[0].ti[0]), &(tr->td[0].count), tr->rdta->numsp, tr->numBranches); 
1320 
1321#ifdef _USE_PTHREADS
1322  masterBarrier(THREAD_EVALUATE_VECTOR, tr);
1323#else
1324  evaluateGenericVectorIterative(tr, 0, tr->cdta->endsite); 
1325#endif
1326}
Note: See TracBrowser for help on using the repository browser.