source: branches/lib/GDE/RAxML/ancestralStates.c

Last change on this file was 13845, checked in by westram, 11 years ago
  • remove executable flag
File size: 22.7 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 <sys/times.h>
33#include <sys/types.h>
34#include <sys/time.h>
35#include <unistd.h>
36#endif
37
38#include <limits.h>
39#include <math.h>
40#include <time.h>
41#include <stdlib.h>
42#include <stdio.h>
43#include <ctype.h>
44#include <string.h>
45
46#include "axml.h"
47
48extern char  workdir[1024];
49extern char run_id[128];
50
51extern const char binaryStateNames[2];   
52extern const char dnaStateNames[4];
53extern const char protStateNames[20];
54extern const char genericStateNames[32];
55
56static char getStateCharacter(int dataType, int state)
57{
58  char 
59    result;
60
61  assert(MIN_MODEL < dataType && dataType < MAX_MODEL);
62
63  switch(dataType)
64    {
65    case BINARY_DATA:
66      result = binaryStateNames[state];
67      break;
68    case DNA_DATA:
69       result = dnaStateNames[state];
70      break;
71    case AA_DATA:
72      result =  protStateNames[state];
73      break; 
74    case GENERIC_32:
75      result = genericStateNames[state];
76      break;
77    default:
78      assert(0);
79    }
80
81  return  result;
82}
83
84static void makeP_Flex_Ancestral(double *rptr, double *EI,  double *EIGN, int numberOfCategories, double *left, const int numStates)
85{
86  int 
87    i,
88    j,
89    k;
90 
91  const int
92    rates = numStates - 1,
93    statesSquare = numStates * numStates;
94
95  double 
96    z1 = 0.0,
97    lz1[64],
98    d1[64];
99
100  assert(numStates <= 64);
101     
102  for(i = 0; i < rates; i++)   
103    lz1[i] = EIGN[i] * z1;
104     
105
106  for(i = 0; i < numberOfCategories; i++)
107    {
108      for(j = 0; j < rates; j++)       
109        d1[j] = EXP (rptr[i] * lz1[j]);
110         
111      for(j = 0; j < numStates; j++)
112        {
113          left[statesSquare * i  + numStates * j] = 1.0;         
114
115          for(k = 1; k < numStates; k++)           
116            left[statesSquare * i + numStates * j + k]  = d1[k-1] * EI[rates * j + (k-1)];           
117        }
118    } 
119}
120
121static void ancestralCat(double *v, double *sumBuffer, double *diagptable, int i, int numStates)
122{
123  int 
124    l,
125    j;
126 
127  double 
128    *ancestral = &sumBuffer[numStates * i],
129    sum = 0.0,
130    *term = (double*)rax_malloc(sizeof(double) * numStates);
131
132  for(l = 0; l < numStates; l++)
133    {
134      double 
135        ump_x1 = 0.0;
136     
137      for(j = 0; j < numStates; j++)   
138        ump_x1 += v[j] * diagptable[l * numStates + j];
139
140      sum += ump_x1;
141      term[l] = ump_x1;
142     
143    }
144               
145  for(l = 0; l < numStates; l++)         
146    ancestral[l] = term[l] / sum;       
147   
148  rax_free(term);
149}
150
151static void newviewFlexCat_Ancestral(int tipCase, double *extEV,
152                                     int *cptr,
153                                     double *x1, double *x2, double *tipVector,
154                                     unsigned char *tipX1, unsigned char *tipX2,
155                                     int n, double *left, double *right, const int numStates, double *diagptable, double *sumBuffer)
156{
157  double   
158    *x3 = (double *)rax_malloc(sizeof(double) * numStates),
159    *le, *ri, *v, *vl, *vr,
160    ump_x1, ump_x2, x1px2;
161  int 
162    i, l, j, scale;
163
164  const int 
165    statesSquare = numStates * numStates;
166
167  switch(tipCase)
168    {
169    case TIP_TIP:
170      {
171        for (i = 0; i < n; i++)
172          {
173            le = &left[cptr[i] * statesSquare];
174            ri = &right[cptr[i] * statesSquare];
175
176            vl = &(tipVector[numStates * tipX1[i]]);
177            vr = &(tipVector[numStates * tipX2[i]]);               
178
179            v  = x3;
180
181            for(l = 0; l < numStates; l++)
182              v[l] = 0.0;
183
184            for(l = 0; l < numStates; l++)
185              {
186                ump_x1 = 0.0;
187                ump_x2 = 0.0;
188
189                for(j = 0; j < numStates; j++)
190                  {
191                    ump_x1 += vl[j] * le[l * numStates + j];
192                    ump_x2 += vr[j] * ri[l * numStates + j];
193                  }
194
195                x1px2 = ump_x1 * ump_x2;
196
197                for(j = 0; j < numStates; j++)           
198                  v[j] += x1px2 * extEV[l * numStates + j];             
199              } 
200                   
201            ancestralCat(v, sumBuffer, &diagptable[cptr[i] * statesSquare], i, numStates);
202          }
203      }
204      break;
205    case TIP_INNER:
206      {
207        for (i = 0; i < n; i++)
208          {
209            le = &left[cptr[i] * statesSquare];
210            ri = &right[cptr[i] * statesSquare];
211
212            vl = &(tipVector[numStates * tipX1[i]]);
213            vr = &x2[numStates * i];
214            v  = x3;
215
216            for(l = 0; l < numStates; l++)
217              v[l] = 0.0;
218
219            for(l = 0; l < numStates; l++)
220              {
221                ump_x1 = 0.0;
222                ump_x2 = 0.0;
223
224                for(j = 0; j < numStates; j++)
225                  {
226                    ump_x1 += vl[j] * le[l * numStates + j];
227                    ump_x2 += vr[j] * ri[l * numStates + j];
228                  }
229
230                x1px2 = ump_x1 * ump_x2;
231
232                for(j = 0; j < numStates; j++)
233                  v[j] += x1px2 * extEV[l * numStates + j];
234              }
235
236            scale = 1;
237            for(l = 0; scale && (l < numStates); l++)
238              scale = ((v[l] < minlikelihood) && (v[l] > minusminlikelihood));     
239
240            if(scale)
241              {
242                for(l = 0; l < numStates; l++)
243                  v[l] *= twotothe256;                   
244              }
245           
246            ancestralCat(v, sumBuffer, &diagptable[cptr[i] * statesSquare], i, numStates);         
247          }
248      }
249      break;
250    case INNER_INNER:
251      for(i = 0; i < n; i++)
252        {
253          le = &left[cptr[i] * statesSquare];
254          ri = &right[cptr[i] * statesSquare];
255
256          vl = &x1[numStates * i];
257          vr = &x2[numStates * i];
258          v = x3;
259
260          for(l = 0; l < numStates; l++)
261            v[l] = 0.0;
262
263          for(l = 0; l < numStates; l++)
264            {
265              ump_x1 = 0.0;
266              ump_x2 = 0.0;
267
268              for(j = 0; j < numStates; j++)
269                {
270                  ump_x1 += vl[j] * le[l * numStates + j];
271                  ump_x2 += vr[j] * ri[l * numStates + j];
272                }
273
274              x1px2 =  ump_x1 * ump_x2;
275
276              for(j = 0; j < numStates; j++)
277                v[j] += x1px2 * extEV[l * numStates + j];
278            }
279
280           scale = 1;
281           for(l = 0; scale && (l < numStates); l++)
282             scale = ((v[l] < minlikelihood) && (v[l] > minusminlikelihood));
283         
284           if(scale)
285             {
286               for(l = 0; l < numStates; l++)
287                 v[l] *= twotothe256;               
288             }
289         
290           ancestralCat(v, sumBuffer, &diagptable[cptr[i] * statesSquare], i, numStates);       
291        }
292      break;
293    default:
294      assert(0);
295    }
296
297  rax_free(x3);
298
299 
300
301}
302
303
304
305static void ancestralGamma(double *_v, double *sumBuffer, double *diagptable, int i, int numStates, int gammaStates)
306{
307  int 
308    statesSquare = numStates * numStates,
309    k,
310    j,
311    l;
312
313  double 
314    *v,
315    *ancestral = &sumBuffer[gammaStates * i],
316    sum = 0.0,
317    *term = (double*)rax_malloc(sizeof(double) * numStates);                 
318
319  for(l = 0; l < numStates; l++)
320    term[l] = 0.0;
321
322  for(k = 0; k < 4; k++)
323    {     
324      v =  &(_v[numStates * k]);
325
326      for(l = 0; l < numStates; l++)
327        {
328          double
329            al = 0.0;
330         
331          for(j = 0; j < numStates; j++)           
332            al += v[j] * diagptable[k * statesSquare + l * numStates + j];
333         
334          term[l] += al;
335          sum += al;
336        }
337    }
338 
339  for(l = 0; l < numStates; l++)       
340    ancestral[l] = term[l] / sum;       
341   
342  rax_free(term);
343}
344
345
346static void newviewFlexGamma_Ancestral(int tipCase,
347                                       double *x1, double *x2, double *extEV, double *tipVector,
348                                       unsigned char *tipX1, unsigned char *tipX2,
349                                       int n, double *left, double *right,
350                                       const int numStates, double *diagptable, double *sumBuffer)
351{
352  double 
353    *v,
354    *x3 = (double*)rax_malloc(sizeof(double) * 4 * numStates);
355  double x1px2;
356  int  i, j, l, k, scale;
357  double *vl, *vr, al, ar;
358
359
360
361  const int 
362    statesSquare = numStates * numStates,
363    gammaStates = 4 * numStates;
364
365  switch(tipCase)
366    {
367    case TIP_TIP:
368      {
369        for(i = 0; i < n; i++)
370          {
371            for(k = 0; k < 4; k++)
372              {
373                vl = &(tipVector[numStates * tipX1[i]]);
374                vr = &(tipVector[numStates * tipX2[i]]);
375                v =  &(x3[numStates * k]);
376
377                for(l = 0; l < numStates; l++)
378                  v[l] = 0;
379
380                for(l = 0; l < numStates; l++)
381                  {
382                    al = 0.0;
383                    ar = 0.0;
384                    for(j = 0; j < numStates; j++)
385                      {
386                        al += vl[j] * left[k * statesSquare + l * numStates + j];
387                        ar += vr[j] * right[k * statesSquare + l * numStates + j];
388                      }
389
390                    x1px2 = al * ar;
391                    for(j = 0; j < numStates; j++)
392                      v[j] += x1px2 * extEV[numStates * l + j];
393                  }     
394              }
395           
396            ancestralGamma(x3, sumBuffer, diagptable, i, numStates, gammaStates);         
397          }
398      }
399      break;
400    case TIP_INNER:
401      {
402        for (i = 0; i < n; i++)
403          {
404            for(k = 0; k < 4; k++)
405              {
406                vl = &(tipVector[numStates * tipX1[i]]);
407                vr = &(x2[gammaStates * i + numStates * k]);
408                v =  &(x3[numStates * k]);
409
410                for(l = 0; l < numStates; l++)
411                  v[l] = 0;
412
413                for(l = 0; l < numStates; l++)
414                  {
415                    al = 0.0;
416                    ar = 0.0;
417                    for(j = 0; j < numStates; j++)
418                      {
419                        al += vl[j] * left[k * statesSquare + l * numStates + j];
420                        ar += vr[j] * right[k * statesSquare + l * numStates + j];
421                      }
422
423                    x1px2 = al * ar;
424                    for(j = 0; j < numStates; j++)
425                      v[j] += x1px2 * extEV[numStates * l + j];
426                  }
427              }
428           
429            v = x3;
430            scale = 1;
431            for(l = 0; scale && (l < gammaStates); l++)
432              scale = (ABS(v[l]) <  minlikelihood);
433
434            if(scale)
435              {
436                for(l = 0; l < gammaStates; l++)
437                  v[l] *= twotothe256;     
438              }
439
440            ancestralGamma(x3, sumBuffer, diagptable, i, numStates, gammaStates);                   
441          }
442      }
443      break;
444    case INNER_INNER:
445      for (i = 0; i < n; i++)
446       {
447         for(k = 0; k < 4; k++)
448           {
449             vl = &(x1[gammaStates * i + numStates * k]);
450             vr = &(x2[gammaStates * i + numStates * k]);
451             v =  &(x3[numStates * k]);
452
453             for(l = 0; l < numStates; l++)
454               v[l] = 0;
455
456             for(l = 0; l < numStates; l++)
457               {
458                 al = 0.0;
459                 ar = 0.0;
460                 for(j = 0; j < numStates; j++)
461                   {
462                     al += vl[j] * left[k * statesSquare + l * numStates + j];
463                     ar += vr[j] * right[k * statesSquare + l * numStates + j];
464                   }
465
466                 x1px2 = al * ar;
467                 for(j = 0; j < numStates; j++)
468                   v[j] += x1px2 * extEV[numStates * l + j];
469               }
470           }
471         
472         v = x3;
473         scale = 1;
474         for(l = 0; scale && (l < gammaStates); l++)
475           scale = ((ABS(v[l]) <  minlikelihood));
476
477         if (scale)
478           {
479             for(l = 0; l < gammaStates; l++)
480               v[l] *= twotothe256;                 
481           }
482         
483         ancestralGamma(x3, sumBuffer, diagptable, i, numStates, gammaStates);           
484       }
485      break;
486    default:
487      assert(0);
488    }
489
490 
491
492  rax_free(x3);
493
494}
495void newviewIterativeAncestral(tree *tr)
496{
497  traversalInfo
498    *ti   = tr->td[0].ti;
499 
500  int 
501    i, 
502    model;
503 
504  assert(!tr->saveMemory);
505 
506  for(i = 1; i < tr->td[0].count; i++)
507    {
508      traversalInfo
509        *tInfo = &ti[i];
510
511      for(model = 0; model < tr->NumberOfModels; model++)
512        {                   
513          double
514            *x1_start = (double*)NULL,
515            *x2_start = (double*)NULL,     
516            *left     = (double*)NULL,
517            *right    = (double*)NULL,                         
518            qz, 
519            rz;
520                                         
521          unsigned char
522            *tipX1 = (unsigned char *)NULL,
523            *tipX2 = (unsigned char *)NULL;
524         
525          size_t           
526            states = (size_t)tr->partitionData[model].states,
527            width = tr->partitionData[model].width;
528                       
529                         
530                 
531          switch(tInfo->tipCase)
532            {
533            case TIP_TIP:                 
534              tipX1    = tr->partitionData[model].yVector[tInfo->qNumber];
535              tipX2    = tr->partitionData[model].yVector[tInfo->rNumber];                                                                                   
536              break;
537            case TIP_INNER:             
538              tipX1    = tr->partitionData[model].yVector[tInfo->qNumber];               
539              x2_start = tr->partitionData[model].xVector[tInfo->rNumber - tr->mxtips - 1];                 
540              break;
541            case INNER_INNER:                                               
542              x1_start = tr->partitionData[model].xVector[tInfo->qNumber - tr->mxtips - 1];
543              x2_start = tr->partitionData[model].xVector[tInfo->rNumber - tr->mxtips - 1];                       
544              break;
545            default:
546              assert(0);
547            }
548         
549          left  = tr->partitionData[model].left;
550          right = tr->partitionData[model].right;
551         
552          if(tr->multiBranch)
553            {
554              qz = tInfo->qz[model];
555              rz = tInfo->rz[model];
556            }
557          else
558            {
559              qz = tInfo->qz[0];
560              rz = tInfo->rz[0];
561            }                                                           
562         
563          switch(tr->rateHetModel)
564            {
565            case CAT:
566              { 
567                double
568                  *diagptable = (double*)rax_malloc(tr->partitionData[model].numberOfCategories * states * states * sizeof(double));
569               
570                makeP_Flex(qz, rz, tr->partitionData[model].perSiteRates,
571                           tr->partitionData[model].EI,
572                           tr->partitionData[model].EIGN,
573                           tr->partitionData[model].numberOfCategories, left, right, states);
574               
575
576                makeP_Flex_Ancestral(tr->partitionData[model].perSiteRates,
577                                     tr->partitionData[model].EI,
578                                     tr->partitionData[model].EIGN,
579                                     tr->partitionData[model].numberOfCategories, diagptable, states);
580                                     
581
582                newviewFlexCat_Ancestral(tInfo->tipCase,  tr->partitionData[model].EV, tr->partitionData[model].rateCategory,
583                                         x1_start, x2_start, tr->partitionData[model].tipVector,
584                                         tipX1, tipX2, width, left, right, states, diagptable, 
585                                         tr->partitionData[model].sumBuffer);
586
587                rax_free(diagptable);
588              }
589              break;
590            case GAMMA:
591            case GAMMA_I:
592              { 
593                double
594                  *diagptable = (double*)rax_malloc(4 * states * states * sizeof(double));
595               
596                makeP_Flex(qz, rz, tr->partitionData[model].gammaRates,
597                           tr->partitionData[model].EI,
598                           tr->partitionData[model].EIGN,
599                           4, left, right, states);
600               
601                makeP_Flex_Ancestral(tr->partitionData[model].gammaRates,
602                                     tr->partitionData[model].EI,
603                                     tr->partitionData[model].EIGN,
604                                     4, diagptable, states);
605               
606
607                newviewFlexGamma_Ancestral(tInfo->tipCase,
608                                           x1_start, x2_start,
609                                           tr->partitionData[model].EV,
610                                           tr->partitionData[model].tipVector,
611                                           tipX1, tipX2,
612                                           width, left, right, states, diagptable, tr->partitionData[model].sumBuffer);
613               
614                rax_free(diagptable);
615              }
616              break;
617            default:
618              assert(0);
619            }                                                                   
620        }   
621    }
622
623}
624static void traversalInfoAncestralRoot(nodeptr p, traversalInfo *ti, int *counter, int maxTips, int numBranches)
625{
626  int 
627    i;
628 
629  nodeptr
630    q = p,
631    r = p->back;
632
633  for(i = 0; i < numBranches; i++)
634    {
635      double 
636        z = sqrt(q->z[i]);
637       if(z < zmin) 
638         z = zmin;
639       if(z > zmax)
640         z = zmax;
641       
642       z = log(z);
643       
644       ti[*counter].qz[i] = z;
645       ti[*counter].rz[i] = z;     
646    }
647
648  if(isTip(r->number, maxTips) && isTip(q->number, maxTips))
649    {
650      ti[*counter].tipCase = TIP_TIP;     
651      ti[*counter].qNumber = q->number;
652      ti[*counter].rNumber = r->number;
653     
654      *counter = *counter + 1;
655    }
656  else
657    {
658      if(isTip(r->number, maxTips) || isTip(q->number, maxTips))
659        {
660          nodeptr tmp;
661         
662          if(isTip(r->number, maxTips))
663            {
664              tmp = r;
665              r = q;
666              q = tmp;
667            }
668         
669          ti[*counter].tipCase = TIP_INNER;       
670          ti[*counter].qNumber = q->number;
671          ti[*counter].rNumber = r->number;               
672         
673          *counter = *counter + 1;
674        }
675      else
676        {         
677          ti[*counter].tipCase = INNER_INNER;       
678          ti[*counter].qNumber = q->number;
679          ti[*counter].rNumber = r->number;             
680         
681          *counter = *counter + 1;
682        }
683    }
684}
685
686
687void newviewGenericAncestral (tree *tr, nodeptr p, boolean atRoot)
688{ 
689  if(atRoot)
690    {     
691      tr->td[0].count = 1;
692      traversalInfoAncestralRoot(p, &(tr->td[0].ti[0]), &(tr->td[0].count), tr->mxtips, tr->numBranches);
693     
694      if(tr->td[0].count > 1)
695        {
696#ifdef _USE_PTHREADS
697          masterBarrier(THREAD_NEWVIEW_ANCESTRAL, tr);
698#else
699          newviewIterativeAncestral(tr);
700#endif
701        }
702    }   
703  else
704    {
705      if(isTip(p->number, tr->mxtips))
706        return;
707     
708     
709      tr->td[0].count = 1;
710      computeTraversalInfo(p, &(tr->td[0].ti[0]), &(tr->td[0].count), tr->mxtips, tr->numBranches);
711     
712      if(tr->td[0].count > 1)
713        {
714#ifdef _USE_PTHREADS
715          masterBarrier(THREAD_NEWVIEW_ANCESTRAL, tr);
716#else
717          newviewIterativeAncestral(tr);
718#endif 
719        }
720    }
721}
722
723typedef struct {
724  double *probs;
725  char c;
726  int states;
727} ancestralState;
728
729
730static void computeAncestralRec(tree *tr, nodeptr p, int *counter, FILE *probsFile, FILE *statesFile, boolean atRoot)
731{
732#ifdef _USE_PTHREADS
733  size_t 
734    accumulatedOffset = 0;
735#endif
736
737  int 
738    model,
739    globalIndex = 0;
740 
741  ancestralState
742    *a = (ancestralState *)rax_malloc(sizeof(ancestralState) * tr->cdta->endsite),
743    *unsortedA = (ancestralState *)rax_malloc(sizeof(ancestralState) * tr->rdta->sites);
744 
745  if(!atRoot)
746    {
747      if(isTip(p->number, tr->mxtips))
748        return;
749 
750      computeAncestralRec(tr, p->next->back,       counter, probsFile, statesFile, atRoot);
751      computeAncestralRec(tr, p->next->next->back, counter, probsFile, statesFile, atRoot);
752
753      newviewGeneric(tr, p);
754    }
755
756  newviewGenericAncestral(tr, p, atRoot);
757
758#ifdef _USE_PTHREADS
759  masterBarrier(THREAD_GATHER_ANCESTRAL, tr);
760#endif
761
762  if(atRoot)
763    {
764      fprintf(probsFile, "ROOT\n");
765      fprintf(statesFile, "ROOT ");
766    }
767  else
768    {
769      fprintf(probsFile, "%d\n", p->number);
770      fprintf(statesFile, "%d ", p->number);
771    }
772
773  for(model = 0; model < tr->NumberOfModels; model++)
774    {
775      int       
776        offset,
777        i,
778        width = tr->partitionData[model].upper - tr->partitionData[model].lower,       
779        states = tr->partitionData[model].states;
780#ifdef _USE_PTHREADS
781      double
782        *ancestral = &tr->ancestralStates[accumulatedOffset];
783#else
784      double 
785        *ancestral = tr->partitionData[model].sumBuffer;
786#endif
787
788      if(tr->rateHetModel == CAT)
789        offset = 1;
790      else
791        offset = 4;           
792
793      for(i = 0; i < width; i++, globalIndex++)
794        {
795          double
796            equal = 1.0 / (double)states,
797            max = -1.0;
798           
799          boolean
800            approximatelyEqual = TRUE;
801
802          int
803            max_l = -1,
804            l;
805         
806          char 
807            c;
808
809          a[globalIndex].states = states;
810          a[globalIndex].probs = (double *)rax_malloc(sizeof(double) * states);
811         
812          for(l = 0; l < states; l++)
813            {
814              double 
815                value = ancestral[offset * states * i + l];
816
817              if(value > max)
818                {
819                  max = value;
820                  max_l = l;
821                }
822             
823              approximatelyEqual = approximatelyEqual && (ABS(equal - value) < 0.000001);
824             
825              a[globalIndex].probs[l] = value;               
826            }
827
828         
829          if(approximatelyEqual)
830            c = '?';     
831          else
832            c = getStateCharacter(tr->partitionData[model].dataType, max_l);
833         
834          a[globalIndex].c = c;   
835        }
836
837#ifdef _USE_PTHREADS
838      accumulatedOffset += width * offset * states;
839#endif           
840    }
841
842  {
843    int 
844      j, 
845      k;
846   
847    for(j = 0; j < tr->cdta->endsite; j++)
848      {
849        for(k = 0; k < tr->rdta->sites; k++)       
850          if(j == tr->patternPosition[k])               
851            {
852              int 
853                sorted = j,
854                unsorted = tr->columnPosition[k] - 1;
855             
856              unsortedA[unsorted].states = a[sorted].states;
857              unsortedA[unsorted].c = a[sorted].c;
858              unsortedA[unsorted].probs = (double*)rax_malloc(sizeof(double) * unsortedA[unsorted].states);
859              memcpy(unsortedA[unsorted].probs,  a[sorted].probs, sizeof(double) * a[sorted].states);         
860            }     
861        } 
862
863    for(k = 0; k < tr->rdta->sites; k++)
864      {
865        for(j = 0; j < unsortedA[k].states; j++)
866          fprintf(probsFile, "%f ", unsortedA[k].probs[j]);
867        fprintf(probsFile, "\n");
868        fprintf(statesFile, "%c", unsortedA[k].c);
869      }
870    fprintf(probsFile, "\n");
871    fprintf(statesFile, "\n");
872  }
873
874
875  *counter = *counter + 1;
876
877  {
878    int j;
879
880    for(j = 0; j < tr->rdta->sites; j++)
881      rax_free(unsortedA[j].probs);
882    for(j = 0; j < tr->cdta->endsite; j++)
883      rax_free(a[j].probs);
884  }
885
886  rax_free(a);
887  rax_free(unsortedA);
888}
889
890static char *ancestralTreeRec(char *treestr, tree *tr, nodeptr p)
891{       
892  if(isTip(p->number, tr->rdta->numsp)) 
893    {         
894      sprintf(treestr, "%s", tr->nameList[p->number]);   
895     
896      while (*treestr) 
897        treestr++;
898    }
899  else 
900    {                   
901      *treestr++ = '(';     
902      treestr = ancestralTreeRec(treestr, tr, p->next->back);     
903      *treestr++ = ',';
904      treestr = ancestralTreeRec(treestr, tr, p->next->next->back);         
905      *treestr++ = ')'; 
906      sprintf(treestr, "%d", p->number);     
907    }
908       
909  while (*treestr) 
910    treestr++;
911 
912  return  treestr;
913}
914
915
916static char *ancestralTree(char *treestr, tree *tr)
917{
918  *treestr++ = '(';
919  treestr = ancestralTreeRec(treestr, tr, tr->leftRootNode);
920  *treestr++ = ',';
921  treestr = ancestralTreeRec(treestr, tr, tr->rightRootNode);   
922  *treestr++ = ')';
923  sprintf(treestr, "ROOT");
924  while (*treestr) 
925    treestr++; 
926 
927  *treestr++ = ';';
928 
929  *treestr++ = '\0';
930  while (*treestr) treestr++;     
931  return  treestr;
932}
933
934void computeAncestralStates(tree *tr, double referenceLikelihood)
935{
936  int 
937    counter = 0;
938 
939  char 
940    treeFileName[2048],
941    ancestralProbsFileName[2048],
942    ancestralStatesFileName[2048];
943
944  FILE
945    *treeFile,
946    *probsFile,
947    *statesFile;
948
949#ifdef _USE_PTHREADS
950  tr->ancestralStates = (double*)rax_malloc(getContiguousVectorLength(tr) * sizeof(double));
951#endif
952
953  strcpy(ancestralProbsFileName,         workdir);
954  strcpy(ancestralStatesFileName,         workdir);
955  strcpy(treeFileName,         workdir);
956
957  strcat(ancestralProbsFileName,         "RAxML_marginalAncestralProbabilities.");
958  strcat(ancestralStatesFileName,        "RAxML_marginalAncestralStates.");
959  strcat(treeFileName,                   "RAxML_nodeLabelledRootedTree.");
960
961  strcat(ancestralProbsFileName,         run_id);
962  strcat(ancestralStatesFileName,        run_id);
963  strcat(treeFileName,                   run_id);
964 
965  probsFile = myfopen(ancestralProbsFileName,   "w");
966  statesFile = myfopen(ancestralStatesFileName, "w");
967  treeFile   = myfopen(treeFileName,            "w");
968
969  assert(tr->leftRootNode == tr->rightRootNode->back);
970 
971  computeAncestralRec(tr, tr->leftRootNode, &counter, probsFile, statesFile, FALSE);
972  computeAncestralRec(tr, tr->rightRootNode, &counter, probsFile, statesFile, FALSE);
973  computeAncestralRec(tr, tr->rightRootNode, &counter, probsFile, statesFile, TRUE);
974 
975  evaluateGeneric(tr, tr->rightRootNode);
976
977  if(fabs(tr->likelihood - referenceLikelihood) > 0.5)
978    {
979      printf("Something suspiciuous is going on with the marginal ancestral probability computations\n");
980      assert(0);
981    } 
982 
983  assert(counter == tr->mxtips - 1);
984   
985  ancestralTree(tr->tree_string, tr);
986
987  fprintf(treeFile, "%s\n", tr->tree_string);
988
989  fclose(probsFile);
990  fclose(statesFile);
991  fclose(treeFile);
992
993  printBothOpen("Marginal Ancestral Probabilities written to file:\n%s\n\n", ancestralProbsFileName);
994  printBothOpen("Ancestral Sequences based on Marginal Ancestral Probabilities written to file:\n%s\n\n", ancestralStatesFileName); 
995  printBothOpen("Node-laballed ROOTED tree written to file:\n%s\n", treeFileName);
996}
997
998
Note: See TracBrowser for help on using the repository browser.