source: trunk/GDE/RAxML/evaluatePartialGenericSpecial.c

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