source: trunk/GDE/RAxML/fastDNAparsimony.c

Last change on this file was 19480, checked in by westram, 17 months ago
File size: 59.5 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
28 *  thousands of taxa and mixed models".
29 *  Bioinformatics 2006; doi: 10.1093/bioinformatics/btl446
30 */
31
32
33#ifndef WIN32
34#include <sys/times.h>
35#include <sys/types.h>
36#include <sys/time.h>
37#include <unistd.h> 
38#endif
39
40#include <limits.h>
41#include <math.h>
42#include <time.h>
43#include <stdlib.h>
44#include <stdio.h>
45#include <ctype.h>
46#include <string.h>
47#include <stdint.h>
48
49#ifdef __AVX
50
51#ifdef __SIM_SSE3
52
53#define _SSE3_WAS_DEFINED
54
55#undef __SIM_SSE3
56
57#endif
58
59#endif
60
61
62#ifdef __SIM_SSE3
63
64#include <xmmintrin.h>
65#include <pmmintrin.h>
66 
67#endif
68
69#ifdef __AVX
70
71#include <xmmintrin.h>
72#include <immintrin.h>
73
74#endif
75
76
77#include "axml.h"
78
79extern const unsigned int mask32[32]; 
80/* vector-specific stuff */
81
82extern char **globalArgv;
83extern int globalArgc;
84
85#ifdef __SIM_SSE3
86
87#define INTS_PER_VECTOR 4
88#define INT_TYPE __m128i
89#define CAST __m128i*
90#define SET_ALL_BITS_ONE _mm_set_epi32(0xFFFFFFFF, 0xFFFFFFFF, 0xFFFFFFFF, 0xFFFFFFFF)
91#define SET_ALL_BITS_ZERO _mm_set_epi32(0x00000000, 0x00000000, 0x00000000, 0x00000000)
92#define VECTOR_LOAD _mm_load_si128
93#define VECTOR_BIT_AND _mm_and_si128
94#define VECTOR_BIT_OR  _mm_or_si128
95#define VECTOR_STORE  _mm_store_si128
96#define VECTOR_AND_NOT _mm_andnot_si128
97
98#endif
99
100#ifdef __AVX
101
102#define INTS_PER_VECTOR 8
103#define INT_TYPE __m256d
104#define CAST double*
105#define SET_ALL_BITS_ONE (__m256d)_mm256_set_epi32(0xFFFFFFFF, 0xFFFFFFFF, 0xFFFFFFFF, 0xFFFFFFFF, 0xFFFFFFFF, 0xFFFFFFFF, 0xFFFFFFFF, 0xFFFFFFFF)
106#define SET_ALL_BITS_ZERO (__m256d)_mm256_set_epi32(0x00000000, 0x00000000, 0x00000000, 0x00000000, 0x00000000, 0x00000000, 0x00000000, 0x00000000)
107#define VECTOR_LOAD _mm256_load_pd
108#define VECTOR_BIT_AND _mm256_and_pd
109#define VECTOR_BIT_OR  _mm256_or_pd
110#define VECTOR_STORE  _mm256_store_pd
111#define VECTOR_AND_NOT _mm256_andnot_pd
112
113#endif
114
115extern double masterTime;
116extern char  workdir[1024];
117extern char run_id[128];
118
119/********************************DNA FUNCTIONS *****************************************************************/
120
121
122
123static void checkSeed(analdef *adef)
124{ 
125  static boolean seedChecked = FALSE;
126
127  if(!seedChecked) 
128    {
129      /*printf("Checking seed\n");*/
130
131      if(adef->parsimonySeed <= 0)
132        {
133          printf("Error: you need to specify a random number seed with \"-p\" for the randomized stepwise addition\n");
134          printf("parsimony algorithm or random tree building algorithm such that runs can be reproduced and debugged ... exiting\n");     
135        }
136 
137      assert(adef->parsimonySeed > 0);
138      seedChecked = TRUE;
139    }
140}
141
142static void getxnodeLocal (nodeptr p)
143{
144  nodeptr  s;
145
146  if((s = p->next)->x || (s = s->next)->x)
147    {
148      p->x = s->x;
149      s->x = 0;
150    }
151}
152
153static void computeTraversalInfoParsimony(nodeptr p, int *ti, int *counter, int maxTips, boolean full)
154{       
155  nodeptr
156    q = p->next->back,
157    r = p->next->next->back;
158 
159  if(! p->x)
160    getxnodeLocal(p); 
161 
162  if(full)
163    {
164       if(q->number > maxTips) 
165         computeTraversalInfoParsimony(q, ti, counter, maxTips, full);
166     
167      if(r->number > maxTips) 
168        computeTraversalInfoParsimony(r, ti, counter, maxTips, full);
169    }
170  else
171    {
172      if(q->number > maxTips && !q->x) 
173        computeTraversalInfoParsimony(q, ti, counter, maxTips, full);
174     
175      if(r->number > maxTips && !r->x) 
176        computeTraversalInfoParsimony(r, ti, counter, maxTips, full);
177    }
178 
179 
180  ti[*counter]     = p->number;
181  ti[*counter + 1] = q->number;
182  ti[*counter + 2] = r->number;
183  *counter = *counter + 4;
184}
185
186
187
188#if (defined(__SIM_SSE3) || defined(__AVX))
189
190static inline unsigned int populationCount(INT_TYPE v_N)
191{
192  unsigned int
193    res[INTS_PER_VECTOR] __attribute__ ((aligned (BYTE_ALIGNMENT)));
194 
195  unsigned int 
196    i,
197    a = 0;
198   
199  VECTOR_STORE((CAST)res, v_N);
200   
201  for(i = 0; i < INTS_PER_VECTOR; i++)
202    a += BIT_COUNT(res[i]);
203   
204  return a;       
205}
206
207#else
208
209static inline unsigned int populationCount(unsigned int n)
210{
211  return BIT_COUNT(n);
212}
213
214#endif
215
216#if (defined(__SIM_SSE3) || defined(__AVX))
217
218void newviewParsimonyIterativeFast(tree *tr)
219{   
220  INT_TYPE
221    allOne = SET_ALL_BITS_ONE;
222
223  int 
224    model,
225    *ti = tr->ti,
226    count = ti[0],
227    index; 
228
229  for(index = 4; index < count; index += 4)
230    {     
231      unsigned int
232        totalScore = 0;
233
234      size_t
235        pNumber = (size_t)ti[index],
236        qNumber = (size_t)ti[index + 1],
237        rNumber = (size_t)ti[index + 2];
238     
239      for(model = 0; model < tr->NumberOfModels; model++)
240        {
241          size_t
242            k,
243            states = tr->partitionData[model].states,
244            width = tr->partitionData[model].parsimonyLength;   
245           
246          unsigned int 
247            i;     
248                 
249          switch(states)
250            {
251            case 2:       
252              {
253                parsimonyNumber
254                  *left[2],
255                  *right[2],
256                  *this[2];
257
258                for(k = 0; k < 2; k++)
259                  {
260                    left[k]  = &(tr->partitionData[model].parsVect[(width * 2 * qNumber) + width * k]);
261                    right[k] = &(tr->partitionData[model].parsVect[(width * 2 * rNumber) + width * k]);
262                    this[k]  = &(tr->partitionData[model].parsVect[(width * 2 * pNumber) + width * k]);
263                  }
264
265                for(i = 0; i < width; i += INTS_PER_VECTOR)
266                  {               
267                    INT_TYPE
268                      s_r, s_l, v_N,
269                      l_A, l_C,
270                      v_A, v_C;         
271                   
272                    s_l = VECTOR_LOAD((CAST)(&left[0][i]));
273                    s_r = VECTOR_LOAD((CAST)(&right[0][i]));
274                    l_A = VECTOR_BIT_AND(s_l, s_r);
275                    v_A = VECTOR_BIT_OR(s_l, s_r);
276                   
277                    s_l = VECTOR_LOAD((CAST)(&left[1][i]));
278                    s_r = VECTOR_LOAD((CAST)(&right[1][i]));
279                    l_C = VECTOR_BIT_AND(s_l, s_r);
280                    v_C = VECTOR_BIT_OR(s_l, s_r);                                                               
281                   
282                    v_N = VECTOR_BIT_OR(l_A, l_C);
283                   
284                    VECTOR_STORE((CAST)(&this[0][i]), VECTOR_BIT_OR(l_A, VECTOR_AND_NOT(v_N, v_A)));
285                    VECTOR_STORE((CAST)(&this[1][i]), VECTOR_BIT_OR(l_C, VECTOR_AND_NOT(v_N, v_C)));                                                                   
286                   
287                    v_N = VECTOR_AND_NOT(v_N, allOne);
288                   
289                    totalScore += populationCount(v_N);           
290                  }
291              }
292              break;
293            case 4:
294              {
295                parsimonyNumber
296                  *left[4],
297                  *right[4],
298                  *this[4];
299
300                for(k = 0; k < 4; k++)
301                  {
302                    left[k]  = &(tr->partitionData[model].parsVect[(width * 4 * qNumber) + width * k]);
303                    right[k] = &(tr->partitionData[model].parsVect[(width * 4 * rNumber) + width * k]);
304                    this[k]  = &(tr->partitionData[model].parsVect[(width * 4 * pNumber) + width * k]);
305                  }
306
307                for(i = 0; i < width; i += INTS_PER_VECTOR)
308                  {               
309                    INT_TYPE
310                      s_r, s_l, v_N,
311                      l_A, l_C, l_G, l_T,
312                      v_A, v_C, v_G, v_T;               
313                   
314                    s_l = VECTOR_LOAD((CAST)(&left[0][i]));
315                    s_r = VECTOR_LOAD((CAST)(&right[0][i]));
316                    l_A = VECTOR_BIT_AND(s_l, s_r);
317                    v_A = VECTOR_BIT_OR(s_l, s_r);
318                   
319                    s_l = VECTOR_LOAD((CAST)(&left[1][i]));
320                    s_r = VECTOR_LOAD((CAST)(&right[1][i]));
321                    l_C = VECTOR_BIT_AND(s_l, s_r);
322                    v_C = VECTOR_BIT_OR(s_l, s_r);
323                   
324                    s_l = VECTOR_LOAD((CAST)(&left[2][i]));
325                    s_r = VECTOR_LOAD((CAST)(&right[2][i]));
326                    l_G = VECTOR_BIT_AND(s_l, s_r);
327                    v_G = VECTOR_BIT_OR(s_l, s_r);
328                   
329                    s_l = VECTOR_LOAD((CAST)(&left[3][i]));
330                    s_r = VECTOR_LOAD((CAST)(&right[3][i]));
331                    l_T = VECTOR_BIT_AND(s_l, s_r);
332                    v_T = VECTOR_BIT_OR(s_l, s_r);
333                   
334                    v_N = VECTOR_BIT_OR(VECTOR_BIT_OR(l_A, l_C), VECTOR_BIT_OR(l_G, l_T));                               
335                   
336                    VECTOR_STORE((CAST)(&this[0][i]), VECTOR_BIT_OR(l_A, VECTOR_AND_NOT(v_N, v_A)));
337                    VECTOR_STORE((CAST)(&this[1][i]), VECTOR_BIT_OR(l_C, VECTOR_AND_NOT(v_N, v_C)));
338                    VECTOR_STORE((CAST)(&this[2][i]), VECTOR_BIT_OR(l_G, VECTOR_AND_NOT(v_N, v_G)));
339                    VECTOR_STORE((CAST)(&this[3][i]), VECTOR_BIT_OR(l_T, VECTOR_AND_NOT(v_N, v_T)));                                                   
340                   
341                    v_N = VECTOR_AND_NOT(v_N, allOne);
342                   
343                    totalScore += populationCount(v_N); 
344                  }
345              }
346              break;
347            case 20:
348              {
349                parsimonyNumber
350                  *left[20],
351                  *right[20],
352                  *this[20];
353
354                for(k = 0; k < 20; k++)
355                  {
356                    left[k]  = &(tr->partitionData[model].parsVect[(width * 20 * qNumber) + width * k]);
357                    right[k] = &(tr->partitionData[model].parsVect[(width * 20 * rNumber) + width * k]);
358                    this[k]  = &(tr->partitionData[model].parsVect[(width * 20 * pNumber) + width * k]);
359                  }
360
361                for(i = 0; i < width; i += INTS_PER_VECTOR)
362                  {               
363                    size_t j;
364                   
365                    INT_TYPE
366                      s_r, s_l, 
367                      v_N = SET_ALL_BITS_ZERO,
368                      l_A[20], 
369                      v_A[20];           
370                   
371                    for(j = 0; j < 20; j++)
372                      {
373                        s_l = VECTOR_LOAD((CAST)(&left[j][i]));
374                        s_r = VECTOR_LOAD((CAST)(&right[j][i]));
375                        l_A[j] = VECTOR_BIT_AND(s_l, s_r);
376                        v_A[j] = VECTOR_BIT_OR(s_l, s_r);
377                       
378                        v_N = VECTOR_BIT_OR(v_N, l_A[j]);
379                      }
380                   
381                    for(j = 0; j < 20; j++)                 
382                      VECTOR_STORE((CAST)(&this[j][i]), VECTOR_BIT_OR(l_A[j], VECTOR_AND_NOT(v_N, v_A[j])));                                                                   
383                   
384                    v_N = VECTOR_AND_NOT(v_N, allOne);
385                   
386                    totalScore += populationCount(v_N);
387                  }
388              }
389              break;
390            default:
391              {
392                parsimonyNumber
393                  *left[32], 
394                  *right[32],
395                  *this[32];
396
397                assert(states <= 32);
398               
399                for(k = 0; k < states; k++)
400                  {
401                    left[k]  = &(tr->partitionData[model].parsVect[(width * states * qNumber) + width * k]);
402                    right[k] = &(tr->partitionData[model].parsVect[(width * states * rNumber) + width * k]);
403                    this[k]  = &(tr->partitionData[model].parsVect[(width * states * pNumber) + width * k]);
404                  }
405
406                for(i = 0; i < width; i += INTS_PER_VECTOR)
407                  {               
408                    size_t j;
409                   
410                    INT_TYPE
411                      s_r, s_l, 
412                      v_N = SET_ALL_BITS_ZERO,
413                      l_A[32], 
414                      v_A[32];           
415                   
416                    for(j = 0; j < states; j++)
417                      {
418                        s_l = VECTOR_LOAD((CAST)(&left[j][i]));
419                        s_r = VECTOR_LOAD((CAST)(&right[j][i]));
420                        l_A[j] = VECTOR_BIT_AND(s_l, s_r);
421                        v_A[j] = VECTOR_BIT_OR(s_l, s_r);
422                       
423                        v_N = VECTOR_BIT_OR(v_N, l_A[j]);
424                      }
425                   
426                    for(j = 0; j < states; j++)             
427                      VECTOR_STORE((CAST)(&this[j][i]), VECTOR_BIT_OR(l_A[j], VECTOR_AND_NOT(v_N, v_A[j])));                                                                   
428                   
429                    v_N = VECTOR_AND_NOT(v_N, allOne);
430                   
431                    totalScore += populationCount(v_N);
432                  }                             
433              }
434            }           
435        }
436
437      tr->parsimonyScore[pNumber] = totalScore + tr->parsimonyScore[rNumber] + tr->parsimonyScore[qNumber];     
438    }
439}
440
441
442
443unsigned int evaluateParsimonyIterativeFast(tree *tr)
444{
445  INT_TYPE
446    allOne = SET_ALL_BITS_ONE;
447
448  size_t 
449    pNumber = (size_t)tr->ti[1],
450    qNumber = (size_t)tr->ti[2];
451
452  int
453    model;
454
455  unsigned int 
456    bestScore = tr->bestParsimony,   
457    sum;
458
459  if(tr->ti[0] > 4)
460    newviewParsimonyIterativeFast(tr); 
461
462  sum = tr->parsimonyScore[pNumber] + tr->parsimonyScore[qNumber];
463
464  for(model = 0; model < tr->NumberOfModels; model++)
465    {
466      size_t
467        k,
468        states = tr->partitionData[model].states,
469        width = tr->partitionData[model].parsimonyLength, 
470        i;
471
472       switch(states)
473         {
474         case 2:
475           {
476             parsimonyNumber
477               *left[2],
478               *right[2];
479             
480             for(k = 0; k < 2; k++)
481               {
482                 left[k]  = &(tr->partitionData[model].parsVect[(width * 2 * qNumber) + width * k]);
483                 right[k] = &(tr->partitionData[model].parsVect[(width * 2 * pNumber) + width * k]);
484               }     
485             
486             for(i = 0; i < width; i += INTS_PER_VECTOR)
487               {                                               
488                 INT_TYPE     
489                   l_A = VECTOR_BIT_AND(VECTOR_LOAD((CAST)(&left[0][i])), VECTOR_LOAD((CAST)(&right[0][i]))),
490                   l_C = VECTOR_BIT_AND(VECTOR_LOAD((CAST)(&left[1][i])), VECTOR_LOAD((CAST)(&right[1][i]))),           
491                   v_N = VECTOR_BIT_OR(l_A, l_C);
492                 
493                 v_N = VECTOR_AND_NOT(v_N, allOne);
494                 
495                 sum += populationCount(v_N);
496                 
497                 if(sum >= bestScore)
498                   return sum;                         
499               }
500           }
501           break;
502         case 4:
503           {
504             parsimonyNumber
505               *left[4],
506               *right[4];
507     
508             for(k = 0; k < 4; k++)
509               {
510                 left[k]  = &(tr->partitionData[model].parsVect[(width * 4 * qNumber) + width * k]);
511                 right[k] = &(tr->partitionData[model].parsVect[(width * 4 * pNumber) + width * k]);
512               }       
513
514             for(i = 0; i < width; i += INTS_PER_VECTOR)
515               {                                               
516                 INT_TYPE     
517                   l_A = VECTOR_BIT_AND(VECTOR_LOAD((CAST)(&left[0][i])), VECTOR_LOAD((CAST)(&right[0][i]))),
518                   l_C = VECTOR_BIT_AND(VECTOR_LOAD((CAST)(&left[1][i])), VECTOR_LOAD((CAST)(&right[1][i]))),
519                   l_G = VECTOR_BIT_AND(VECTOR_LOAD((CAST)(&left[2][i])), VECTOR_LOAD((CAST)(&right[2][i]))),
520                   l_T = VECTOR_BIT_AND(VECTOR_LOAD((CAST)(&left[3][i])), VECTOR_LOAD((CAST)(&right[3][i]))),
521                   v_N = VECTOR_BIT_OR(VECTOR_BIT_OR(l_A, l_C), VECTOR_BIT_OR(l_G, l_T));     
522                 
523                 v_N = VECTOR_AND_NOT(v_N, allOne);
524                 
525                 sum += populationCount(v_N);
526                 
527                 if(sum >= bestScore)           
528                   return sum;         
529               }                 
530           }
531           break;
532         case 20:
533           {
534             parsimonyNumber
535               *left[20],
536               *right[20];
537             
538              for(k = 0; k < 20; k++)
539                {
540                  left[k]  = &(tr->partitionData[model].parsVect[(width * 20 * qNumber) + width * k]);
541                  right[k] = &(tr->partitionData[model].parsVect[(width * 20 * pNumber) + width * k]);
542                } 
543           
544              for(i = 0; i < width; i += INTS_PER_VECTOR)
545                {                             
546                  int 
547                    j;
548                 
549                  INT_TYPE     
550                    l_A,
551                    v_N = SET_ALL_BITS_ZERO;     
552                 
553                  for(j = 0; j < 20; j++)
554                    {
555                      l_A = VECTOR_BIT_AND(VECTOR_LOAD((CAST)(&left[j][i])), VECTOR_LOAD((CAST)(&right[j][i])));
556                      v_N = VECTOR_BIT_OR(l_A, v_N);
557                    }
558                 
559                  v_N = VECTOR_AND_NOT(v_N, allOne);
560                 
561                  sum += populationCount(v_N);         
562                 
563                  if(sum >= bestScore)     
564                    return sum;                       
565                }
566           }
567           break;
568         default:
569           {
570             parsimonyNumber
571               *left[32], 
572               *right[32]; 
573
574             assert(states <= 32);
575
576             for(k = 0; k < states; k++)
577               {
578                 left[k]  = &(tr->partitionData[model].parsVect[(width * states * qNumber) + width * k]);
579                 right[k] = &(tr->partitionData[model].parsVect[(width * states * pNumber) + width * k]);
580               } 
581           
582             for(i = 0; i < width; i += INTS_PER_VECTOR)
583               {                               
584                 size_t
585                   j;
586                 
587                 INT_TYPE     
588                   l_A,
589                   v_N = SET_ALL_BITS_ZERO;     
590                 
591                 for(j = 0; j < states; j++)
592                   {
593                     l_A = VECTOR_BIT_AND(VECTOR_LOAD((CAST)(&left[j][i])), VECTOR_LOAD((CAST)(&right[j][i])));
594                     v_N = VECTOR_BIT_OR(l_A, v_N);
595                   }
596                 
597                 v_N = VECTOR_AND_NOT(v_N, allOne);
598                 
599                 sum += populationCount(v_N);         
600                 
601                 if(sum >= bestScore)         
602                   return sum;                 
603               }
604           }
605         }
606    }
607 
608  return sum;
609}
610
611
612#else
613
614void newviewParsimonyIterativeFast(tree *tr)
615{   
616  int 
617    model,
618    *ti = tr->ti,
619    count = ti[0],
620    index; 
621
622  for(index = 4; index < count; index += 4)
623    {     
624      unsigned int
625        totalScore = 0;
626
627      size_t
628        pNumber = (size_t)ti[index],
629        qNumber = (size_t)ti[index + 1],
630        rNumber = (size_t)ti[index + 2];
631     
632      for(model = 0; model < tr->NumberOfModels; model++)
633        {
634          size_t
635            // k,
636            states = tr->partitionData[model].states,
637            width = tr->partitionData[model].parsimonyLength;   
638           
639          unsigned int 
640            i;     
641                 
642          switch(states)
643            {
644            case 2:       
645              {
646                parsimonyNumber
647                  *left[2],
648                  *right[2],
649                  *this[2];
650               
651                parsimonyNumber
652                   o_A,
653                   o_C,
654                   t_A,
655                   t_C, 
656                   t_N;
657               
658                for(size_t k = 0; k < 2; k++)
659                  {
660                    left[k]  = &(tr->partitionData[model].parsVect[(width * 2 * qNumber) + width * k]);
661                    right[k] = &(tr->partitionData[model].parsVect[(width * 2 * rNumber) + width * k]);
662                    this[k]  = &(tr->partitionData[model].parsVect[(width * 2 * pNumber) + width * k]);
663                  }
664
665                for(i = 0; i < width; i++)
666                  {               
667                    t_A = left[0][i] & right[0][i];
668                    t_C = left[1][i] & right[1][i];               
669
670                    o_A = left[0][i] | right[0][i];
671                    o_C = left[1][i] | right[1][i];
672                 
673                    t_N = ~(t_A | t_C);   
674
675                    this[0][i] = t_A | (t_N & o_A);
676                    this[1][i] = t_C | (t_N & o_C);               
677                   
678                    totalScore += populationCount(t_N);   
679                  }
680              }
681              break;
682            case 4:
683              {
684                parsimonyNumber
685                  *left[4],
686                  *right[4],
687                  *this[4];
688
689                for(size_t k = 0; k < 4; k++)
690                  {
691                    left[k]  = &(tr->partitionData[model].parsVect[(width * 4 * qNumber) + width * k]);
692                    right[k] = &(tr->partitionData[model].parsVect[(width * 4 * rNumber) + width * k]);
693                    this[k]  = &(tr->partitionData[model].parsVect[(width * 4 * pNumber) + width * k]);
694                  }
695
696                parsimonyNumber
697                   o_A,
698                   o_C,
699                   o_G,
700                   o_T,
701                   t_A,
702                   t_C,
703                   t_G,
704                   t_T, 
705                   t_N;
706
707                for(i = 0; i < width; i++)
708                  {               
709                    t_A = left[0][i] & right[0][i];
710                    t_C = left[1][i] & right[1][i];
711                    t_G = left[2][i] & right[2][i];       
712                    t_T = left[3][i] & right[3][i];
713
714                    o_A = left[0][i] | right[0][i];
715                    o_C = left[1][i] | right[1][i];
716                    o_G = left[2][i] | right[2][i];       
717                    o_T = left[3][i] | right[3][i];
718
719                    t_N = ~(t_A | t_C | t_G | t_T);       
720
721                    this[0][i] = t_A | (t_N & o_A);
722                    this[1][i] = t_C | (t_N & o_C);
723                    this[2][i] = t_G | (t_N & o_G);
724                    this[3][i] = t_T | (t_N & o_T); 
725                   
726                    totalScore += populationCount(t_N);   
727                  }
728              }
729              break;
730            case 20:
731              {
732                parsimonyNumber
733                  *left[20],
734                  *right[20],
735                  *this[20];
736
737                parsimonyNumber
738                  o_A[20],
739                  t_A[20],       
740                  t_N;
741
742                for(size_t k = 0; k < 20; k++)
743                  {
744                    left[k]  = &(tr->partitionData[model].parsVect[(width * 20 * qNumber) + width * k]);
745                    right[k] = &(tr->partitionData[model].parsVect[(width * 20 * rNumber) + width * k]);
746                    this[k]  = &(tr->partitionData[model].parsVect[(width * 20 * pNumber) + width * k]);
747                  }
748
749                for(i = 0; i < width; i++)
750                  {               
751                    size_t k;
752                   
753                    t_N = 0;
754
755                    for(k = 0; k < 20; k++)
756                      {
757                        t_A[k] = left[k][i] & right[k][i];
758                        o_A[k] = left[k][i] | right[k][i];
759                        t_N = t_N | t_A[k];
760                      }
761                   
762                    t_N = ~t_N;
763
764                    for(k = 0; k < 20; k++)                   
765                      this[k][i] = t_A[k] | (t_N & o_A[k]);               
766                   
767                    totalScore += populationCount(t_N); 
768                  }
769              }
770              break;
771            default:
772              {         
773                parsimonyNumber
774                  *left[32],
775                  *right[32],
776                  *this[32];
777               
778                parsimonyNumber
779                  o_A[32],
780                  t_A[32],       
781                  t_N;
782               
783                assert(states <= 32);
784               
785                for(size_t k = 0; k < states; k++)
786                  {
787                    left[k]  = &(tr->partitionData[model].parsVect[(width * states * qNumber) + width * k]);
788                    right[k] = &(tr->partitionData[model].parsVect[(width * states * rNumber) + width * k]);
789                    this[k]  = &(tr->partitionData[model].parsVect[(width * states * pNumber) + width * k]);
790                  }
791               
792                for(i = 0; i < width; i++)
793                  {               
794                    t_N = 0;
795                   
796                    for(size_t k = 0; k < states; k++)
797                      {
798                        t_A[k] = left[k][i] & right[k][i];
799                        o_A[k] = left[k][i] | right[k][i];
800                        t_N = t_N | t_A[k];
801                      }
802                   
803                    t_N = ~t_N;
804                   
805                    for(size_t k = 0; k < states; k++)
806                      this[k][i] = t_A[k] | (t_N & o_A[k]);               
807                   
808                    totalScore += populationCount(t_N); 
809                  }
810              }                       
811            } 
812        }
813
814      tr->parsimonyScore[pNumber] = totalScore + tr->parsimonyScore[rNumber] + tr->parsimonyScore[qNumber];     
815    }
816}
817
818
819
820unsigned int evaluateParsimonyIterativeFast(tree *tr)
821{
822  size_t 
823    pNumber = (size_t)tr->ti[1],
824    qNumber = (size_t)tr->ti[2];
825
826  int
827    model;
828
829  unsigned int 
830    bestScore = tr->bestParsimony,   
831    sum;
832
833  if(tr->ti[0] > 4)
834    newviewParsimonyIterativeFast(tr); 
835
836  sum = tr->parsimonyScore[pNumber] + tr->parsimonyScore[qNumber];
837
838  for(model = 0; model < tr->NumberOfModels; model++)
839    {
840      size_t
841        k,
842        states = tr->partitionData[model].states,
843        width = tr->partitionData[model].parsimonyLength, 
844        i;
845
846       switch(states)
847         {
848         case 2:
849           {
850             parsimonyNumber
851               t_A,
852               t_C,           
853               t_N,
854               *left[2],
855               *right[2];
856             
857             for(k = 0; k < 2; k++)
858               {
859                 left[k]  = &(tr->partitionData[model].parsVect[(width * 2 * qNumber) + width * k]);
860                 right[k] = &(tr->partitionData[model].parsVect[(width * 2 * pNumber) + width * k]);
861               }     
862             
863             for(i = 0; i < width; i++)
864               {                                               
865                 t_A = left[0][i] & right[0][i];
866                 t_C = left[1][i] & right[1][i];
867                 
868                  t_N = ~(t_A | t_C);
869
870                  sum += populationCount(t_N);   
871                 
872                 if(sum >= bestScore)
873                   return sum;                         
874               }
875           }
876           break;
877         case 4:
878           {
879             parsimonyNumber
880               t_A,
881               t_C,
882               t_G,
883               t_T,
884               t_N,
885               *left[4],
886               *right[4];
887     
888             for(k = 0; k < 4; k++)
889               {
890                 left[k]  = &(tr->partitionData[model].parsVect[(width * 4 * qNumber) + width * k]);
891                 right[k] = &(tr->partitionData[model].parsVect[(width * 4 * pNumber) + width * k]);
892               }       
893
894             for(i = 0; i < width; i++)
895               {                                               
896                  t_A = left[0][i] & right[0][i];
897                  t_C = left[1][i] & right[1][i];
898                  t_G = left[2][i] & right[2][i];         
899                  t_T = left[3][i] & right[3][i];
900
901                  t_N = ~(t_A | t_C | t_G | t_T);
902
903                  sum += populationCount(t_N);     
904                 
905                 if(sum >= bestScore)           
906                   return sum;         
907               }                 
908           }
909           break;
910         case 20:
911           {
912             parsimonyNumber
913               t_A,
914               t_N,
915               *left[20],
916               *right[20];
917             
918              for(k = 0; k < 20; k++)
919                {
920                  left[k]  = &(tr->partitionData[model].parsVect[(width * 20 * qNumber) + width * k]);
921                  right[k] = &(tr->partitionData[model].parsVect[(width * 20 * pNumber) + width * k]);
922                } 
923           
924              for(i = 0; i < width; i++)
925                { 
926                  t_N = 0;
927                 
928                  for(k = 0; k < 20; k++)
929                    {
930                      t_A = left[k][i] & right[k][i];
931                      t_N = t_N | t_A;
932                    }
933               
934                  t_N = ~t_N;
935
936                  sum += populationCount(t_N);     
937                 
938                  if(sum >= bestScore)     
939                    return sum;                       
940                }
941           }
942           break;
943         default:
944           {
945             parsimonyNumber
946               t_A,
947               t_N,
948               *left[32], 
949               *right[32]; 
950
951             assert(states <= 32);
952
953             for(k = 0; k < states; k++)
954               {
955                 left[k]  = &(tr->partitionData[model].parsVect[(width * states * qNumber) + width * k]);
956                 right[k] = &(tr->partitionData[model].parsVect[(width * states * pNumber) + width * k]);
957               } 
958           
959             for(i = 0; i < width; i++)
960               {                               
961                 t_N = 0;
962                 
963                 for(k = 0; k < states; k++)
964                   {
965                     t_A = left[k][i] & right[k][i];
966                     t_N = t_N | t_A;
967                   }
968               
969                  t_N = ~t_N;
970
971                  sum += populationCount(t_N);     
972                                                 
973                 if(sum >= bestScore)                     
974                   return sum;                     
975               }                     
976           }
977         }
978    }
979 
980  return sum;
981}
982
983#endif
984
985
986
987
988
989
990static unsigned int evaluateParsimony(tree *tr, nodeptr p, boolean full)
991{
992  volatile unsigned int result;
993  nodeptr q = p->back;
994  int
995    *ti = tr->ti,
996    counter = 4;
997 
998  ti[1] = p->number;
999  ti[2] = q->number;
1000
1001  if(full)
1002    {
1003      if(p->number > tr->mxtips)
1004        computeTraversalInfoParsimony(p, ti, &counter, tr->mxtips, full);
1005      if(q->number > tr->mxtips)
1006        computeTraversalInfoParsimony(q, ti, &counter, tr->mxtips, full); 
1007    }
1008  else
1009    {
1010      if(p->number > tr->mxtips && !p->x)
1011        computeTraversalInfoParsimony(p, ti, &counter, tr->mxtips, full);
1012      if(q->number > tr->mxtips && !q->x)
1013        computeTraversalInfoParsimony(q, ti, &counter, tr->mxtips, full); 
1014    }
1015
1016  ti[0] = counter;
1017
1018  result = evaluateParsimonyIterativeFast(tr);
1019
1020  return result;
1021}
1022
1023
1024static void newviewParsimony(tree *tr, nodeptr  p)
1025{     
1026  if(p->number <= tr->mxtips)
1027    return;
1028
1029  {
1030    int 
1031      counter = 4;     
1032           
1033    computeTraversalInfoParsimony(p, tr->ti, &counter, tr->mxtips, FALSE);             
1034    tr->ti[0] = counter;           
1035   
1036    newviewParsimonyIterativeFast(tr);     
1037  }
1038}
1039
1040
1041
1042
1043
1044/****************************************************************************************************************************************/
1045
1046static void insertParsimony (tree *tr, nodeptr p, nodeptr q)
1047{
1048  nodeptr  r;
1049 
1050  r = q->back;
1051 
1052  hookupDefault(p->next,       q, tr->numBranches);
1053  hookupDefault(p->next->next, r, tr->numBranches); 
1054   
1055  newviewParsimony(tr, p);     
1056} 
1057
1058/*
1059  static nodeptr buildNewTip (tree *tr, nodeptr p)
1060  {
1061  nodeptr  q;
1062 
1063  q = tr->nodep[(tr->nextnode)++];
1064  hookupDefault(p, q, tr->numBranches);
1065  q->next->back = (nodeptr)NULL;
1066  q->next->next->back = (nodeptr)NULL;
1067  assert(q == q->next->next->next);
1068  assert(q->x || q->next->x || q->next->next->x);
1069  return  q;
1070  }
1071*/
1072
1073static nodeptr buildNewTip (tree *tr, nodeptr p)
1074{ 
1075  nodeptr  q;
1076
1077  q = tr->nodep[(tr->nextnode)++];
1078  hookupDefault(p, q, tr->numBranches);
1079  q->next->back = (nodeptr)NULL;
1080  q->next->next->back = (nodeptr)NULL;
1081 
1082  return  q;
1083} 
1084
1085static void buildSimpleTree (tree *tr, int ip, int iq, int ir)
1086{   
1087  nodeptr  p, s;
1088  int  i;
1089 
1090  i = MIN(ip, iq);
1091  if (ir < i)  i = ir; 
1092  tr->start = tr->nodep[i];
1093  tr->ntips = 3;
1094  p = tr->nodep[ip];
1095  hookupDefault(p, tr->nodep[iq], tr->numBranches);
1096  s = buildNewTip(tr, tr->nodep[ir]);
1097  insertParsimony(tr, s, p);
1098}
1099
1100
1101static void testInsertParsimony (tree *tr, nodeptr p, nodeptr q)
1102{ 
1103  unsigned int 
1104    mp;
1105 
1106  nodeptr 
1107    r = q->back;   
1108
1109  boolean
1110    doIt = TRUE;
1111   
1112  if(tr->grouped)
1113    {
1114      int 
1115        rNumber = tr->constraintVector[r->number],
1116        qNumber = tr->constraintVector[q->number],
1117        pNumber = tr->constraintVector[p->number];
1118
1119      doIt = FALSE;
1120     
1121      if(pNumber == -9)
1122        pNumber = checker(tr, p->back);
1123      if(pNumber == -9)
1124        doIt = TRUE;
1125      else
1126        {
1127          if(qNumber == -9)
1128            qNumber = checker(tr, q);
1129
1130          if(rNumber == -9)
1131            rNumber = checker(tr, r);
1132
1133          if(pNumber == rNumber || pNumber == qNumber)
1134            doIt = TRUE;       
1135        }
1136    }
1137
1138  if(doIt)
1139    {
1140      insertParsimony(tr, p, q);   
1141 
1142      mp = evaluateParsimony(tr, p->next->next, FALSE);         
1143     
1144      if(mp < tr->bestParsimony)
1145        {
1146          tr->bestParsimony = mp;
1147          tr->insertNode = q;
1148          tr->removeNode = p;
1149        }
1150 
1151      hookupDefault(q, r, tr->numBranches);
1152      p->next->next->back = p->next->back = (nodeptr) NULL;
1153    }
1154       
1155  return;
1156} 
1157
1158
1159static void restoreTreeParsimony(tree *tr, nodeptr p, nodeptr q)
1160{ 
1161  nodeptr
1162    r = q->back;
1163 
1164  int counter = 4;
1165 
1166  hookupDefault(p->next,       q, tr->numBranches);
1167  hookupDefault(p->next->next, r, tr->numBranches);
1168 
1169  computeTraversalInfoParsimony(p, tr->ti, &counter, tr->mxtips, FALSE);             
1170  tr->ti[0] = counter;
1171   
1172  newviewParsimonyIterativeFast(tr); 
1173}
1174
1175
1176static void addTraverseParsimony (tree *tr, nodeptr p, nodeptr q, int mintrav, int maxtrav, boolean doAll)
1177{       
1178  if (doAll || (--mintrav <= 0))               
1179    testInsertParsimony(tr, p, q);                       
1180
1181  if (((q->number > tr->mxtips)) && ((--maxtrav > 0) || doAll))
1182    {         
1183      addTraverseParsimony(tr, p, q->next->back, mintrav, maxtrav, doAll);           
1184      addTraverseParsimony(tr, p, q->next->next->back, mintrav, maxtrav, doAll);                     
1185    }
1186}
1187
1188
1189static nodeptr findAnyTipFast(nodeptr p, int numsp)
1190{ 
1191  return  (p->number <= numsp)? p : findAnyTipFast(p->next->back, numsp);
1192} 
1193
1194
1195static void makePermutationFast(int *perm, int n, analdef *adef)
1196{   
1197  int 
1198    i, 
1199    j, 
1200    k;
1201
1202  checkSeed(adef);     
1203
1204  for (i = 1; i <= n; i++)   
1205    perm[i] = i;               
1206
1207  for (i = 1; i <= n; i++) 
1208    {     
1209      double d =  randum(&adef->parsimonySeed);
1210
1211      k =  (int)((double)(n + 1 - i) * d);
1212     
1213      j        = perm[i];
1214
1215      perm[i]     = perm[i + k];
1216      perm[i + k] = j; 
1217    }
1218}
1219
1220static nodeptr  removeNodeParsimony (nodeptr p, tree *tr)
1221{ 
1222  nodeptr  q, r;         
1223
1224  q = p->next->back;
1225  r = p->next->next->back;   
1226   
1227  hookupDefault(q, r, tr->numBranches);
1228
1229  p->next->next->back = p->next->back = (node *) NULL;
1230 
1231  return  q;
1232}
1233
1234static int rearrangeParsimony(tree *tr, nodeptr p, int mintrav, int maxtrav, boolean doAll) 
1235{   
1236  nodeptr 
1237    p1, 
1238    p2, 
1239    q, 
1240    q1, 
1241    q2;
1242 
1243  int     
1244    mintrav2; 
1245
1246  boolean
1247    doP = TRUE,
1248    doQ = TRUE;
1249           
1250  if (maxtrav > tr->ntips - 3) 
1251    maxtrav = tr->ntips - 3; 
1252
1253  assert(mintrav == 1);
1254
1255  if(maxtrav < mintrav)
1256    return 0;
1257
1258  q = p->back;
1259
1260  if(tr->constrained)
1261    {   
1262      if(! tipHomogeneityChecker(tr, p->back, 0))
1263        doP = FALSE;
1264       
1265      if(! tipHomogeneityChecker(tr, q->back, 0))
1266        doQ = FALSE;
1267                       
1268      if(doQ == FALSE && doP == FALSE)
1269        return 0;
1270    } 
1271
1272  if((p->number > tr->mxtips) && doP) 
1273    {     
1274      p1 = p->next->back;
1275      p2 = p->next->next->back;
1276     
1277      if ((p1->number > tr->mxtips) || (p2->number > tr->mxtips)) 
1278        {                 
1279          removeNodeParsimony(p, tr);           
1280
1281          if ((p1->number > tr->mxtips)) 
1282            {
1283              addTraverseParsimony(tr, p, p1->next->back, mintrav, maxtrav, doAll);         
1284              addTraverseParsimony(tr, p, p1->next->next->back, mintrav, maxtrav, doAll);         
1285            }
1286         
1287          if ((p2->number > tr->mxtips)) 
1288            {
1289              addTraverseParsimony(tr, p, p2->next->back, mintrav, maxtrav, doAll);
1290              addTraverseParsimony(tr, p, p2->next->next->back, mintrav, maxtrav, doAll);         
1291            }
1292           
1293           
1294          hookupDefault(p->next,       p1, tr->numBranches); 
1295          hookupDefault(p->next->next, p2, tr->numBranches);                       
1296
1297          newviewParsimony(tr, p);
1298        }
1299    } 
1300       
1301  if ((q->number > tr->mxtips) && (maxtrav > 0) && doQ) 
1302    {
1303      q1 = q->next->back;
1304      q2 = q->next->next->back;
1305
1306      if (
1307          (
1308           (q1->number > tr->mxtips) && 
1309           ((q1->next->back->number > tr->mxtips) || (q1->next->next->back->number > tr->mxtips))
1310           )
1311          ||
1312          (
1313           (q2->number > tr->mxtips) && 
1314           ((q2->next->back->number > tr->mxtips) || (q2->next->next->back->number > tr->mxtips))
1315           )
1316          )
1317        {         
1318
1319          removeNodeParsimony(q, tr);
1320         
1321          mintrav2 = mintrav > 2 ? mintrav : 2;
1322         
1323          if ((q1->number > tr->mxtips)) 
1324            {
1325              addTraverseParsimony(tr, q, q1->next->back, mintrav2 , maxtrav, doAll);
1326              addTraverseParsimony(tr, q, q1->next->next->back, mintrav2 , maxtrav, doAll);         
1327            }
1328         
1329          if ((q2->number > tr->mxtips)) 
1330            {
1331              addTraverseParsimony(tr, q, q2->next->back, mintrav2 , maxtrav, doAll);
1332              addTraverseParsimony(tr, q, q2->next->next->back, mintrav2 , maxtrav, doAll);         
1333            }     
1334           
1335          hookupDefault(q->next,       q1, tr->numBranches); 
1336          hookupDefault(q->next->next, q2, tr->numBranches);
1337           
1338          newviewParsimony(tr, q);
1339        }
1340    }
1341
1342  return 1;
1343} 
1344
1345
1346static void restoreTreeRearrangeParsimony(tree *tr)
1347{   
1348  removeNodeParsimony(tr->removeNode, tr); 
1349  restoreTreeParsimony(tr, tr->removeNode, tr->insertNode); 
1350}
1351
1352/*
1353static boolean isInformative2(tree *tr, int site)
1354{
1355  int
1356    informativeCounter = 0,
1357    check[256],   
1358    j,   
1359    undetermined = 15;
1360
1361  unsigned char
1362    nucleotide,
1363    target = 0;
1364       
1365  for(j = 0; j < 256; j++)
1366    check[j] = 0;
1367 
1368  for(j = 1; j <= tr->mxtips; j++)
1369    {     
1370      nucleotide = tr->yVector[j][site];           
1371      check[nucleotide] =  check[nucleotide] + 1;                 
1372    }
1373 
1374 
1375  if(check[1] > 1)
1376    {
1377      informativeCounter++;   
1378      target = target | 1;
1379    }
1380  if(check[2] > 1)
1381    {
1382      informativeCounter++;
1383      target = target | 2;
1384    }
1385  if(check[4] > 1)
1386    {
1387      informativeCounter++;
1388      target = target | 4;
1389    }
1390  if(check[8] > 1)
1391    {
1392      informativeCounter++;
1393      target = target | 8;
1394    }
1395         
1396  if(informativeCounter >= 2)
1397    return TRUE;   
1398  else
1399    {       
1400      for(j = 0; j < undetermined; j++)
1401        {
1402          if(j == 3 || j == 5 || j == 6 || j == 7 || j == 9 || j == 10 || j == 11 ||
1403             j == 12 || j == 13 || j == 14)
1404            {
1405              if(check[j] > 1)
1406                {
1407                  if(!(target & j))
1408                    return TRUE;
1409                }
1410            }
1411        }
1412    }
1413     
1414  return FALSE;     
1415}
1416*/
1417
1418static boolean isInformative(tree *tr, int dataType, int site)
1419{
1420  int
1421    informativeCounter = 0,
1422    check[256],   
1423    j,   
1424    undetermined = getUndetermined(dataType);
1425
1426  const unsigned int
1427    *bitVector = getBitVector(dataType);
1428
1429  unsigned char
1430    nucleotide;
1431 
1432       
1433  for(j = 0; j < 256; j++)
1434    check[j] = 0;
1435 
1436  for(j = 1; j <= tr->mxtips; j++)
1437    {     
1438      nucleotide = tr->yVector[j][site];           
1439      check[nucleotide] =  check[nucleotide] + 1;
1440      assert(bitVector[nucleotide] > 0);                   
1441    }
1442 
1443  for(j = 0; j < undetermined; j++)
1444    {
1445      if(check[j] > 0)
1446        informativeCounter++;   
1447    } 
1448         
1449  if(informativeCounter <= 1)
1450    return FALSE;   
1451  else
1452    {       
1453      for(j = 0; j < undetermined; j++)
1454        {
1455          if(check[j] > 1)
1456            return TRUE;
1457        } 
1458    }
1459     
1460  return FALSE;     
1461}
1462
1463
1464static void determineUninformativeSites(tree *tr, int *informative)
1465{
1466  int 
1467    i,
1468    number = 0;
1469
1470  /*
1471     Not all characters are useful in constructing a parsimony tree.
1472     Invariant characters, those that have the same state in all taxa,
1473     are obviously useless and are ignored by the method. Characters in
1474     which a state occurs in only one taxon are also ignored.
1475     All these characters are called parsimony uninformative.
1476
1477     Alternative definition: informative columns contain at least two types
1478     of nucleotides, and each nucleotide must appear at least twice in each
1479     column. Kind of a pain if we intend to check for this when using, e.g.,
1480     amibiguous DNA encoding.
1481  */
1482
1483  for(i = 0; i < tr->cdta->endsite; i++)
1484    {
1485      if(isInformative(tr, tr->dataVector[i], i))
1486        informative[i] = 1;
1487      else
1488        {
1489          informative[i] = 0;
1490          number++;
1491        }           
1492    }
1493 
1494 
1495  /* printf("Uninformative Patterns: %d\n", number); */
1496}
1497
1498
1499static void reorderNodes(tree *tr, nodeptr *np, nodeptr p, int *count)
1500{
1501  int i, found = 0;
1502
1503  if((p->number <= tr->mxtips))   
1504    return;
1505  else
1506    {             
1507      for(i = tr->mxtips + 1; (i <= (tr->mxtips + tr->mxtips - 1)) && (found == 0); i++)
1508        {
1509          if (p == np[i] || p == np[i]->next || p == np[i]->next->next)
1510            {
1511              if(p == np[i])                           
1512                tr->nodep[*count + tr->mxtips + 1] = np[i];                             
1513              else
1514                {
1515                  if(p == np[i]->next)           
1516                    tr->nodep[*count + tr->mxtips + 1] = np[i]->next;                     
1517                  else             
1518                    tr->nodep[*count + tr->mxtips + 1] = np[i]->next->next;                                 
1519                }
1520
1521              found = 1;                     
1522              *count = *count + 1;
1523            }
1524        }           
1525     
1526      assert(found != 0);
1527
1528      reorderNodes(tr, np, p->next->back, count);     
1529      reorderNodes(tr, np, p->next->next->back, count);               
1530    }
1531}
1532
1533
1534
1535
1536
1537 
1538static void compressDNA(tree *tr, int *informative, boolean saveMemory)
1539{
1540  size_t
1541    totalNodes,
1542    i,
1543    model;
1544 
1545  if(saveMemory)
1546    totalNodes = (size_t)tr->innerNodes + 1 + (size_t)tr->mxtips;
1547  else
1548    totalNodes = 2 * (size_t)tr->mxtips;
1549
1550 
1551
1552  for(model = 0; model < (size_t) tr->NumberOfModels; model++)
1553    {
1554      size_t
1555        k,
1556        states = (size_t)tr->partitionData[model].states,       
1557        compressedEntries,
1558        compressedEntriesPadded,
1559        entries = 0, 
1560        lower = tr->partitionData[model].lower,
1561        upper = tr->partitionData[model].upper;
1562
1563      parsimonyNumber
1564        **compressedTips = (parsimonyNumber **)rax_malloc(states * sizeof(parsimonyNumber*)),
1565        *compressedValues = (parsimonyNumber *)rax_malloc(states * sizeof(parsimonyNumber));
1566     
1567      for(i = lower; i < upper; i++)   
1568        if(informative[i])
1569          entries += (size_t)tr->cdta->aliaswgt[i];     
1570 
1571      compressedEntries = entries / PCF;
1572
1573      if(entries % PCF != 0)
1574        compressedEntries++;
1575
1576#if (defined(__SIM_SSE3) || defined(__AVX))
1577      if(compressedEntries % INTS_PER_VECTOR != 0)
1578        compressedEntriesPadded = compressedEntries + (INTS_PER_VECTOR - (compressedEntries % INTS_PER_VECTOR));
1579      else
1580        compressedEntriesPadded = compressedEntries;
1581#else
1582      compressedEntriesPadded = compressedEntries;
1583#endif     
1584
1585     
1586      tr->partitionData[model].parsVect = (parsimonyNumber *)rax_malloc((size_t)compressedEntriesPadded * states * totalNodes * sizeof(parsimonyNumber));
1587     
1588      for(i = 0; i < compressedEntriesPadded * states * totalNodes; i++)     
1589        tr->partitionData[model].parsVect[i] = 0;         
1590
1591      for(i = 0; i < (size_t)tr->mxtips; i++)
1592        {
1593          size_t
1594            w = 0,
1595            compressedIndex = 0,
1596            compressedCounter = 0,
1597            index = 0;
1598
1599          for(k = 0; k < states; k++)
1600            {
1601              compressedTips[k] = &(tr->partitionData[model].parsVect[(compressedEntriesPadded * states * (i + 1)) + (compressedEntriesPadded * k)]);
1602              compressedValues[k] = 0;
1603            }               
1604             
1605          for(index = lower; index < (size_t)upper; index++)
1606            {
1607              if(informative[index])
1608                {
1609                  const unsigned int 
1610                    *bitValue = getBitVector(tr->partitionData[model].dataType);
1611
1612                  parsimonyNumber
1613                    value = bitValue[tr->yVector[i + 1][index]];         
1614             
1615                  for(w = 0; w < (size_t)tr->cdta->aliaswgt[index]; w++)
1616                    {     
1617                      for(k = 0; k < states; k++)
1618                        {
1619                          if(value & mask32[k])
1620                            compressedValues[k] |= mask32[compressedCounter];
1621                        }
1622                     
1623                      compressedCounter++;
1624                 
1625                      if(compressedCounter == PCF)
1626                        {
1627                          for(k = 0; k < states; k++)
1628                            {
1629                              compressedTips[k][compressedIndex] = compressedValues[k];
1630                              compressedValues[k] = 0;
1631                            }                   
1632                         
1633                          compressedCounter = 0;
1634                          compressedIndex++;
1635                        }
1636                    }
1637                }
1638            }
1639                           
1640          for(;compressedIndex < compressedEntriesPadded; compressedIndex++)
1641            {   
1642              for(;compressedCounter < PCF; compressedCounter++)             
1643                for(k = 0; k < states; k++)
1644                  compressedValues[k] |= mask32[compressedCounter];               
1645         
1646              for(k = 0; k < states; k++)
1647                {
1648                  compressedTips[k][compressedIndex] = compressedValues[k];
1649                  compressedValues[k] = 0;
1650                }                     
1651             
1652              compressedCounter = 0;
1653            }           
1654        }               
1655 
1656      tr->partitionData[model].parsimonyLength = compressedEntriesPadded;   
1657
1658      rax_free(compressedTips);
1659      rax_free(compressedValues);
1660    }
1661 
1662  tr->parsimonyScore = (unsigned int*)rax_malloc(sizeof(unsigned int) * totalNodes); 
1663         
1664  for(i = 0; i < totalNodes; i++) 
1665    tr->parsimonyScore[i] = 0;
1666}
1667
1668
1669
1670static void stepwiseAddition(tree *tr, nodeptr p, nodeptr q)
1671{           
1672  nodeptr
1673    r = q->back;
1674
1675  unsigned int 
1676    mp;
1677 
1678  int 
1679    counter = 4;
1680 
1681  p->next->back = q;
1682  q->back = p->next;
1683
1684  p->next->next->back = r;
1685  r->back = p->next->next;
1686   
1687  computeTraversalInfoParsimony(p, tr->ti, &counter, tr->mxtips, FALSE);             
1688  tr->ti[0] = counter;
1689  tr->ti[1] = p->number;
1690  tr->ti[2] = p->back->number;
1691   
1692  mp = evaluateParsimonyIterativeFast(tr);
1693 
1694  if(mp < tr->bestParsimony)
1695    {   
1696      tr->bestParsimony = mp;
1697      tr->insertNode = q;     
1698    }
1699 
1700  q->back = r;
1701  r->back = q;
1702   
1703  if(q->number > tr->mxtips && tr->parsimonyScore[q->number] > 0)
1704    {         
1705      stepwiseAddition(tr, p, q->next->back);         
1706      stepwiseAddition(tr, p, q->next->next->back);                         
1707    }
1708}
1709
1710static void markNodesInTree(nodeptr p, tree *tr, unsigned char *nodesInTree)
1711{
1712  if(isTip(p->number, tr->mxtips))
1713    nodesInTree[p->number] = 1;
1714  else
1715    {
1716      markNodesInTree(p->next->back, tr, nodesInTree);
1717      markNodesInTree(p->next->next->back, tr, nodesInTree);
1718    }
1719
1720}
1721
1722void makeParsimonyTreeFast(tree *tr, analdef *adef, boolean full)
1723{   
1724  nodeptr 
1725    p, 
1726    f;   
1727
1728  size_t
1729    model;
1730
1731  int 
1732    i, 
1733    nextsp,
1734    *perm        = (int *)rax_malloc((size_t)(tr->mxtips + 1) * sizeof(int)),
1735    *informative = (int *)rax_malloc(sizeof(int) * (size_t)tr->cdta->endsite); 
1736
1737  unsigned int 
1738    randomMP, 
1739    startMP;       
1740
1741  /* double t; */
1742
1743  determineUninformativeSites(tr, informative);     
1744
1745  compressDNA(tr, informative, FALSE);
1746
1747  rax_free(informative); 
1748
1749  tr->ti = (int*)rax_malloc(sizeof(int) * 4 * (size_t)tr->mxtips); 
1750 
1751  /*t = gettime();*/
1752
1753  if(!full)
1754    {                   
1755      unsigned char
1756        *nodesInTree = (unsigned char*)rax_calloc((size_t)(tr->mxtips + 1), sizeof(unsigned char));           
1757       
1758      tr->start = findAnyTipFast(tr->start, tr->rdta->numsp);
1759       
1760      tr->bestParsimony = INT_MAX;
1761
1762      evaluateParsimony(tr, tr->start->back, TRUE);
1763               
1764      assert(tr->start);
1765     
1766      checkSeed(adef);
1767
1768      markNodesInTree(tr->start, tr, nodesInTree);
1769      markNodesInTree(tr->start->back, tr, nodesInTree);
1770       
1771      int j1 = tr->ntips + 1;
1772       
1773      if(tr->grouped)
1774        {
1775          for(i = 1; i <= tr->mxtips; i++)     
1776            {
1777              if(tr->constraintVector[i] == -1) 
1778                {
1779                  perm[j1++] = i;
1780                  tr->constraintVector[i] = -9;
1781                }
1782            }
1783        }
1784      else
1785        {
1786          if(tr->constrained)
1787            { 
1788              for(i = 1; i <= tr->mxtips; i++)
1789                tr->constraintVector[i] = 0;
1790             
1791              for(i = 1; i <= tr->mxtips; i++)
1792                {                 
1793                  if(nodesInTree[i] == 0)             
1794                    perm[j1++] = i;
1795                  else
1796                    tr->constraintVector[i] = 1;                   
1797                }
1798            }
1799          else
1800            {
1801              for(i = 1; i <= tr->mxtips; i++)     
1802                if(nodesInTree[i] == 0) 
1803                  perm[j1++] = i;
1804            }
1805        }
1806       
1807      for(i = tr->ntips + 1; i <= tr->mxtips; i++) 
1808        {           
1809          int k, j;
1810         
1811          k =  (int)((double)(tr->mxtips + 1 - i) * randum(&adef->parsimonySeed));
1812         
1813          assert(i + k <= tr->mxtips);
1814          j        = perm[i];
1815          perm[i]     = perm[i + k];
1816          perm[i + k] = j;
1817        }   
1818       
1819      f = tr->start;     
1820
1821      rax_free(nodesInTree);
1822    }
1823  else
1824    {
1825      assert(!tr->constrained);
1826
1827      makePermutationFast(perm, tr->mxtips, adef);
1828     
1829      tr->ntips = 0;   
1830     
1831      tr->nextnode = tr->mxtips + 1;       
1832     
1833      buildSimpleTree(tr, perm[1], perm[2], perm[3]);     
1834     
1835      f = tr->start;
1836    }     
1837 
1838  while(tr->ntips < tr->mxtips) 
1839    {   
1840      nodeptr q;
1841     
1842      tr->bestParsimony = INT_MAX;
1843      nextsp = ++(tr->ntips);             
1844      p = tr->nodep[perm[nextsp]];                 
1845      q = tr->nodep[(tr->nextnode)++];
1846      p->back = q;
1847      q->back = p;
1848       
1849      if(tr->grouped && !full)
1850        {
1851          int 
1852            number = p->back->number;           
1853
1854          tr->constraintVector[number] = -9;
1855        }
1856         
1857      stepwiseAddition(tr, q, f->back);                 
1858     
1859      {
1860        nodeptr   
1861          r = tr->insertNode->back;
1862       
1863        int counter = 4;
1864       
1865        hookupDefault(q->next,       tr->insertNode, tr->numBranches);
1866        hookupDefault(q->next->next, r, tr->numBranches);
1867       
1868        computeTraversalInfoParsimony(q, tr->ti, &counter, tr->mxtips, FALSE);             
1869        tr->ti[0] = counter;
1870       
1871        newviewParsimonyIterativeFast(tr);     
1872      }
1873    }   
1874 
1875  //printf("ADD: %d\n", tr->bestParsimony);
1876 
1877  nodeRectifier(tr);
1878 
1879  if(adef->stepwiseAdditionOnly == FALSE)
1880    {
1881      randomMP = tr->bestParsimony;       
1882     
1883      do
1884        {
1885          startMP = randomMP;
1886          nodeRectifier(tr);
1887          for(i = 1; i <= tr->mxtips + tr->mxtips - 2; i++)
1888            {
1889              rearrangeParsimony(tr, tr->nodep[i], 1, 20, FALSE);
1890              if(tr->bestParsimony < randomMP)
1891                {               
1892                  restoreTreeRearrangeParsimony(tr);
1893                  randomMP = tr->bestParsimony;
1894                }
1895            }                             
1896        }
1897      while(randomMP < startMP);
1898    }
1899 
1900  //printf("OPT: %d\n", tr->bestParsimony);
1901
1902 
1903     
1904  rax_free(perm); 
1905  rax_free(tr->parsimonyScore);
1906 
1907  for(model = 0; model < (size_t) tr->NumberOfModels; model++)
1908    rax_free(tr->partitionData[model].parsVect);
1909 
1910  rax_free(tr->ti);
1911} 
1912
1913
1914
1915 
1916static void insertRandom (nodeptr p, nodeptr q, int numBranches)
1917{
1918  nodeptr  r;
1919 
1920  r = q->back;
1921 
1922  hookupDefault(p->next,       q, numBranches);
1923  hookupDefault(p->next->next, r, numBranches); 
1924} 
1925
1926
1927
1928
1929
1930
1931
1932static void buildSimpleTreeRandom (tree *tr, int ip, int iq, int ir)
1933{   
1934  nodeptr  p, s;
1935  int  i;
1936 
1937  i = MIN(ip, iq);
1938  if (ir < i)  i = ir; 
1939  tr->start = tr->nodep[i];
1940  tr->ntips = 3;
1941  p = tr->nodep[ip];
1942  hookupDefault(p, tr->nodep[iq], tr->numBranches);
1943  s = buildNewTip(tr, tr->nodep[ir]);
1944  insertRandom(s, p, tr->numBranches);
1945}
1946
1947int checker(tree *tr, nodeptr p)
1948{
1949  int group = tr->constraintVector[p->number];
1950
1951  if(isTip(p->number, tr->mxtips))
1952    {
1953      group = tr->constraintVector[p->number];
1954      return group;
1955    }
1956  else
1957    {
1958      if(group != -9) 
1959        return group;
1960
1961      group = checker(tr, p->next->back);
1962      if(group != -9) 
1963        return group;
1964
1965      group = checker(tr, p->next->next->back);
1966      if(group != -9) 
1967        return group;
1968
1969      return -9;
1970    }
1971}
1972
1973
1974
1975
1976
1977
1978
1979static int markBranches(nodeptr *branches, nodeptr p, int *counter, int numsp)
1980{
1981  if(isTip(p->number, numsp))
1982    return 0;
1983  else
1984    {
1985      branches[*counter] = p->next;
1986      branches[*counter + 1] = p->next->next;
1987     
1988      *counter = *counter + 2;
1989     
1990      return ((2 + markBranches(branches, p->next->back, counter, numsp) + 
1991               markBranches(branches, p->next->next->back, counter, numsp)));
1992    }
1993}
1994
1995
1996
1997nodeptr findAnyTip(nodeptr p, int numsp)
1998{ 
1999  return  isTip(p->number, numsp) ? p : findAnyTip(p->next->back, numsp);
2000} 
2001
2002
2003int randomInt(int n)
2004{
2005  return rand() %n;
2006}
2007
2008void makePermutation(int *perm, int n, analdef *adef)
2009{   
2010  int  i, j, k;
2011
2012  checkSeed(adef);         
2013
2014  for (i = 1; i <= n; i++)   
2015    perm[i] = i;               
2016
2017  for (i = 1; i <= n; i++) 
2018    {   
2019      k =  (int)((double)(n + 1 - i) * randum(&adef->parsimonySeed));
2020
2021      assert(i + k <= n);
2022     
2023      j        = perm[i];
2024      perm[i]     = perm[i + k];
2025      perm[i + k] = j; 
2026    }
2027}
2028
2029
2030
2031
2032
2033
2034
2035
2036boolean tipHomogeneityChecker(tree *tr, nodeptr p, int grouping)
2037{
2038  if(isTip(p->number, tr->mxtips))
2039    {
2040      if(tr->constraintVector[p->number] != grouping) 
2041        return FALSE;
2042      else 
2043        return TRUE;
2044    }
2045  else
2046    {   
2047      return  (tipHomogeneityChecker(tr, p->next->back, grouping) && tipHomogeneityChecker(tr, p->next->next->back,grouping));     
2048    }
2049}
2050
2051
2052
2053
2054
2055
2056
2057void makeRandomTree(tree *tr, analdef *adef)
2058{ 
2059  nodeptr p, f, randomBranch;   
2060  int nextsp;
2061  int *perm, branchCounter;
2062  nodeptr *branches;
2063 
2064  branches = (nodeptr *)rax_malloc(sizeof(nodeptr) * (2 * tr->mxtips));
2065  perm = (int *)rax_malloc((tr->mxtips + 1) * sizeof(int));                         
2066 
2067  makePermutation(perm, tr->mxtips, adef);             
2068 
2069  tr->ntips = 0;               
2070  tr->nextnode = tr->mxtips + 1;   
2071 
2072  buildSimpleTreeRandom(tr, perm[1], perm[2], perm[3]);
2073 
2074  while (tr->ntips < tr->mxtips) 
2075    {         
2076      tr->bestParsimony = INT_MAX;
2077      nextsp = ++(tr->ntips);             
2078      p = tr->nodep[perm[nextsp]];
2079     
2080      /*printf("ADDING SPECIES %d\n", nextsp);*/
2081     
2082      buildNewTip(tr, p);       
2083     
2084      f = findAnyTip(tr->start, tr->mxtips);
2085      f = f->back;
2086     
2087      branchCounter = 1;
2088      branches[0] = f;
2089      markBranches(branches, f, &branchCounter, tr->mxtips);
2090
2091      assert(branchCounter == ((2 * (tr->ntips - 1)) - 3));
2092     
2093      randomBranch = branches[randomInt(branchCounter)];
2094     
2095      insertRandom(p->back, randomBranch, tr->numBranches);
2096     
2097    }
2098  rax_free(perm);           
2099  rax_free(branches);
2100}
2101
2102
2103
2104void nodeRectifier(tree *tr)
2105{
2106  nodeptr *np = (nodeptr *)rax_malloc(2 * tr->mxtips * sizeof(nodeptr));
2107  int i;
2108  int count = 0;
2109 
2110  tr->start       = tr->nodep[1];
2111  tr->rooted      = FALSE;
2112
2113  /* TODO why is tr->rooted set to FALSE here ?*/
2114 
2115  for(i = tr->mxtips + 1; i <= (tr->mxtips + tr->mxtips - 1); i++)
2116    np[i] = tr->nodep[i];           
2117 
2118  reorderNodes(tr, np, tr->start->back, &count); 
2119
2120 
2121  rax_free(np);
2122}
2123
2124
2125void makeParsimonyTree(tree *tr, analdef *adef)
2126{ 
2127  makeParsimonyTreeFast(tr, adef, TRUE);       
2128}
2129
2130void makeParsimonyTreeIncomplete(tree *tr, analdef *adef)
2131{   
2132  makeParsimonyTreeFast(tr, adef, FALSE);       
2133}
2134 
2135
2136static void setupBranchMetaInfo(tree *tr, nodeptr p, int nTips, branchInfo *bInf)
2137{
2138  int 
2139    countBranches = tr->branchCounter;
2140
2141  if(isTip(p->number, tr->mxtips))   
2142    {     
2143      p->bInf       = &bInf[countBranches];
2144      p->back->bInf = &bInf[countBranches];                           
2145
2146      bInf[countBranches].oP = p;
2147      bInf[countBranches].oQ = p->back;
2148     
2149      bInf[countBranches].epa->leftNodeNumber = p->number;
2150      bInf[countBranches].epa->rightNodeNumber = p->back->number;
2151         
2152      bInf[countBranches].epa->branchNumber = countBranches;                     
2153      bInf[countBranches].epa->originalBranchLength = p->z[0];
2154
2155      tr->branchCounter =  tr->branchCounter + 1;
2156      return;
2157    }
2158  else
2159    {
2160      nodeptr q;
2161      assert(p == p->next->next->next);
2162
2163      p->bInf       = &bInf[countBranches];
2164      p->back->bInf = &bInf[countBranches];
2165
2166      bInf[countBranches].oP = p;
2167      bInf[countBranches].oQ = p->back;
2168
2169      bInf[countBranches].epa->leftNodeNumber = p->number;
2170      bInf[countBranches].epa->rightNodeNumber = p->back->number;
2171
2172               
2173      bInf[countBranches].epa->branchNumber = countBranches;
2174      bInf[countBranches].epa->originalBranchLength = p->z[0];
2175
2176      tr->branchCounter =  tr->branchCounter + 1;     
2177
2178      q = p->next;
2179
2180      while(q != p)
2181        {
2182          setupBranchMetaInfo(tr, q->back, nTips, bInf);       
2183          q = q->next;
2184        }
2185     
2186      return;
2187    }
2188}
2189 
2190
2191
2192static void setupJointFormat(tree *tr, nodeptr p, int ntips, branchInfo *bInf, int *count)
2193{
2194  if(isTip(p->number, tr->mxtips))   
2195    {     
2196      p->bInf->epa->jointLabel = *count;
2197      *count = *count + 1;
2198           
2199      return;
2200    }
2201  else
2202    {                           
2203      setupJointFormat(tr, p->next->back, ntips, bInf, count);           
2204      setupJointFormat(tr, p->next->next->back, ntips, bInf, count);     
2205     
2206      p->bInf->epa->jointLabel = *count;
2207      *count = *count + 1; 
2208     
2209      return;
2210    }
2211}
2212 
2213
2214
2215
2216
2217
2218static void setupBranchInfo(tree *tr, nodeptr q)
2219{
2220  nodeptr
2221    originalNode = tr->nodep[tr->mxtips + 1];
2222
2223  int 
2224    count = 0;
2225
2226  tr->branchCounter = 0;
2227
2228  setupBranchMetaInfo(tr, q, tr->ntips, tr->bInf);
2229   
2230  assert(tr->branchCounter == tr->numberOfBranches);
2231
2232  if(tr->wasRooted)
2233    {
2234      assert(tr->leftRootNode->back == tr->rightRootNode);
2235      assert(tr->leftRootNode       == tr->rightRootNode->back);     
2236
2237      if(!isTip(tr->leftRootNode->number, tr->mxtips))
2238        {
2239          setupJointFormat(tr,  tr->leftRootNode->next->back, tr->ntips, tr->bInf, &count);
2240          setupJointFormat(tr,  tr->leftRootNode->next->next->back, tr->ntips, tr->bInf, &count);
2241        }
2242     
2243       tr->leftRootNode->bInf->epa->jointLabel = count;
2244       tr->rootLabel = count;
2245       count = count + 1;
2246
2247       if(!isTip(tr->rightRootNode->number, tr->mxtips))
2248         {
2249          setupJointFormat(tr,  tr->rightRootNode->next->back, tr->ntips, tr->bInf, &count);
2250          setupJointFormat(tr,  tr->rightRootNode->next->next->back, tr->ntips, tr->bInf, &count);
2251        }             
2252    }
2253  else
2254    {
2255      setupJointFormat(tr, originalNode->back, tr->ntips, tr->bInf, &count);
2256      setupJointFormat(tr, originalNode->next->back, tr->ntips, tr->bInf, &count);
2257      setupJointFormat(tr, originalNode->next->next->back, tr->ntips, tr->bInf, &count);     
2258    } 
2259
2260  assert(count == tr->numberOfBranches);
2261}
2262
2263static void testInsertFast(tree *tr, nodeptr r, nodeptr q)
2264{
2265  unsigned int
2266    result;
2267 
2268  nodeptr 
2269    x = q->back;     
2270 
2271  int 
2272    i,
2273    *inserts = tr->inserts;
2274                   
2275  assert(!tr->grouped);                             
2276 
2277  hookupDefault(r->next,       q, tr->numBranches);
2278  hookupDefault(r->next->next, x, tr->numBranches);                             
2279   
2280  newviewParsimony(tr, r);   
2281   
2282  for(i = 0; i < tr->numberOfTipsForInsertion; i++)
2283    {                               
2284      hookupDefault(r, tr->nodep[inserts[i]], tr->numBranches);
2285     
2286      tr->bestParsimony = INT_MAX;
2287
2288      result = evaluateParsimony(tr, r, FALSE);           
2289     
2290      r->back = (nodeptr) NULL;
2291      tr->nodep[inserts[i]]->back = (nodeptr) NULL;
2292     
2293      tr->bInf[q->bInf->epa->branchNumber].epa->parsimonyScore[i] = result;                     
2294    }
2295 
2296  hookupDefault(q, x, tr->numBranches);
2297 
2298  r->next->next->back = r->next->back = (nodeptr) NULL;
2299 
2300}
2301
2302
2303static void traverseTree(tree *tr, nodeptr r, nodeptr q)
2304{       
2305  testInsertFast(tr, r, q);
2306
2307  if(!isTip(q->number, tr->rdta->numsp))
2308    {   
2309      nodeptr
2310        a = q->next;
2311
2312      while(a != q)
2313        {
2314          traverseTree(tr, r, a->back);
2315          a = a->next;
2316        }     
2317    }
2318} 
2319
2320
2321
2322typedef struct
2323  {
2324    unsigned int parsimonyScore; 
2325    int number;
2326  }
2327  infoMP;
2328
2329
2330static int infoCompare(const void *p1, const void *p2)
2331{
2332  infoMP *rc1 = (infoMP *)p1;
2333  infoMP *rc2 = (infoMP *)p2;
2334
2335  unsigned int i = rc1->parsimonyScore;
2336  unsigned int j = rc2->parsimonyScore;
2337
2338  if (i > j)
2339    return (1);
2340  if (i < j)
2341    return (-1);
2342  return (0);
2343}
2344
2345void classifyMP(tree *tr, analdef *adef)
2346{   
2347  int 
2348    *informative = (int *)rax_malloc(sizeof(int) * (size_t)tr->cdta->endsite),
2349    j, 
2350    *perm;   
2351 
2352  infoMP
2353    *inf;
2354   
2355  nodeptr     
2356    r, 
2357    q;   
2358
2359  char   
2360    jointFormatTreeFileName[1024], 
2361    originalLabelledTreeFileName[1024],
2362    labelledTreeFileName[1024],
2363    likelihoodWeightsFileName[1024];
2364             
2365  FILE   
2366    *likelihoodWeightsFile,
2367    *treeFile;
2368 
2369  unsigned int 
2370    score;       
2371
2372  assert(adef->restart);
2373
2374  determineUninformativeSites(tr, informative);     
2375
2376  compressDNA(tr, informative, TRUE);
2377
2378  rax_free(informative); 
2379
2380  tr->ti = (int*)rax_malloc(sizeof(int) * 4 * (size_t)tr->mxtips);   
2381 
2382  tr->numberOfBranches = 2 * tr->ntips - 3;
2383
2384  printBothOpen("\nRAxML Evolutionary Placement Algorithm using parsimony\n"); 
2385
2386  tr->bestParsimony = INT_MAX;
2387
2388  score = evaluateParsimony(tr, tr->start->back, TRUE);
2389
2390  printBothOpen("\nparsimony score of reference tree: %u\n\n", score);
2391
2392  perm        = (int *)rax_calloc(((size_t)tr->mxtips) + 1, sizeof(int));
2393  tr->inserts = (int *)rax_calloc((size_t)tr->mxtips,   sizeof(int));
2394
2395  markTips(tr->start,       perm, tr->mxtips);
2396  markTips(tr->start->back, perm ,tr->mxtips);
2397
2398  tr->numberOfTipsForInsertion = 0;
2399
2400  for(int i = 1; i <= tr->mxtips; i++)
2401    {
2402      if(perm[i] == 0)
2403        {
2404          tr->inserts[tr->numberOfTipsForInsertion] = i;
2405          tr->numberOfTipsForInsertion = tr->numberOfTipsForInsertion + 1;
2406        }
2407    }   
2408
2409  rax_free(perm);
2410 
2411  printBothOpen("RAxML will place %d Query Sequences into the %d branches of the reference tree with %d taxa\n\n",  tr->numberOfTipsForInsertion, (2 * tr->ntips - 3), tr->ntips); 
2412
2413  assert(tr->numberOfTipsForInsertion == (tr->mxtips - tr->ntips));     
2414
2415  tr->bInf              = (branchInfo*)rax_malloc(tr->numberOfBranches * sizeof(branchInfo)); 
2416
2417  for(int i = 0; i < tr->numberOfBranches; i++)
2418    {     
2419      tr->bInf[i].epa = (epaBranchData*)rax_malloc(sizeof(epaBranchData));                   
2420
2421      tr->bInf[i].epa->parsimonyScore = (unsigned int*)rax_malloc(tr->numberOfTipsForInsertion * sizeof(unsigned int));
2422
2423      for(j = 0; j < tr->numberOfTipsForInsertion; j++)
2424        tr->bInf[i].epa->parsimonyScore[j] = INT_MAX;
2425               
2426      tr->bInf[i].epa->branchNumber = i;
2427      tr->bInf[i].epa->countThem = (int*)rax_calloc(tr->numberOfTipsForInsertion, sizeof(int));
2428     
2429      sprintf(tr->bInf[i].epa->branchLabel, "I%d", i);     
2430    } 
2431
2432  r = tr->nodep[(tr->nextnode)++]; 
2433   
2434
2435  q = findAnyTip(tr->start, tr->rdta->numsp);
2436
2437  assert(isTip(q->number, tr->rdta->numsp));
2438  assert(!isTip(q->back->number, tr->rdta->numsp));
2439         
2440  q = q->back; 
2441 
2442  setupBranchInfo(tr, q);   
2443 
2444  traverseTree(tr, r, q);
2445                   
2446  printBothOpen("Overall Classification time: %f\n\n", gettime() - masterTime); 
2447 
2448 
2449  strcpy(jointFormatTreeFileName,      workdir); 
2450  strcpy(originalLabelledTreeFileName, workdir); 
2451  strcpy(labelledTreeFileName,         workdir); 
2452  strcpy(likelihoodWeightsFileName,    workdir);
2453
2454  strcat(jointFormatTreeFileName,      "RAxML_portableTree."); 
2455  strcat(originalLabelledTreeFileName, "RAxML_originalLabelledTree.");
2456  strcat(labelledTreeFileName,         "RAxML_labelledTree.");
2457  strcat(likelihoodWeightsFileName,    "RAxML_equallyParsimoniousPlacements.");
2458 
2459  strcat(jointFormatTreeFileName,      run_id); 
2460  strcat(originalLabelledTreeFileName, run_id);
2461  strcat(labelledTreeFileName,         run_id);
2462  strcat(likelihoodWeightsFileName,    run_id);
2463
2464  strcat(jointFormatTreeFileName,      ".jplace");
2465 
2466  rax_free(tr->tree_string);
2467
2468  tr->treeStringLength *= 16;
2469
2470  tr->tree_string  = (char*)rax_calloc(tr->treeStringLength, sizeof(char)); 
2471 
2472 
2473
2474  treeFile = myfopen(originalLabelledTreeFileName, "wb");
2475  Tree2StringClassify(tr->tree_string, tr, tr->inserts, TRUE, FALSE, FALSE);
2476  fprintf(treeFile, "%s\n", tr->tree_string);   
2477  fclose(treeFile); 
2478
2479  treeFile = myfopen(jointFormatTreeFileName, "wb");
2480  Tree2StringClassify(tr->tree_string, tr, tr->inserts, TRUE, TRUE, FALSE);
2481 
2482  fprintf(treeFile, "{\n");
2483  fprintf(treeFile, "\t\"tree\": \"%s\", \n", tr->tree_string);
2484  fprintf(treeFile, "\t\"placements\": [\n");
2485     
2486
2487  inf = (infoMP*)rax_malloc(sizeof(infoMP) * tr->numberOfBranches); 
2488 
2489  /* joint format */
2490       
2491 
2492
2493  for(int i = 0; i < tr->numberOfTipsForInsertion; i++)
2494    {
2495      unsigned int
2496        lmax;
2497     
2498      int         
2499        validEntries = tr->numberOfBranches;
2500     
2501      for(j =  0; j < tr->numberOfBranches; j++) 
2502        {
2503          inf[j].parsimonyScore = tr->bInf[j].epa->parsimonyScore[i];     
2504          inf[j].number         = tr->bInf[j].epa->jointLabel;
2505        }
2506     
2507      qsort(inf, tr->numberOfBranches, sizeof(infoMP), infoCompare);               
2508     
2509      j = 0;
2510     
2511      lmax = inf[0].parsimonyScore;
2512     
2513      fprintf(treeFile, "\t{\"p\":[");
2514     
2515      while(j < validEntries && inf[j].parsimonyScore == lmax)   
2516        {           
2517          if(j > 0)
2518            {
2519              if(tr->wasRooted && inf[j].number == tr->rootLabel)                 
2520                assert(0);               
2521              else
2522                fprintf(treeFile, ",[%d, %u]", inf[j].number, inf[j].parsimonyScore);
2523            }
2524          else
2525            {
2526              if(tr->wasRooted && inf[j].number == tr->rootLabel)                 
2527                assert(0);               
2528              else
2529                fprintf(treeFile, "[%d, %u]", inf[j].number, inf[j].parsimonyScore);
2530            }     
2531         
2532          j++;
2533        }
2534     
2535      if(i == tr->numberOfTipsForInsertion - 1)
2536        fprintf(treeFile, "], \"n\":[\"%s\"]}\n", tr->nameList[tr->inserts[i]]);
2537      else
2538        fprintf(treeFile, "], \"n\":[\"%s\"]},\n", tr->nameList[tr->inserts[i]]);
2539    } 
2540   
2541  fprintf(treeFile, "\t ],\n");
2542  /*  fprintf(treeFile, "\t\"metadata\": {\"invocation\": \"RAxML EPA parsimony\"},\n");*/
2543  fprintf(treeFile, "\t\"metadata\": {\"invocation\": ");
2544
2545  fprintf(treeFile, "\"");
2546 
2547  {
2548    int i;
2549   
2550    for(i = 0; i < globalArgc; i++)
2551      fprintf(treeFile,"%s ", globalArgv[i]);
2552  }
2553  fprintf(treeFile, "\", \"raxml_version\": \"%s\"", programVersion);
2554  fprintf(treeFile,"},\n");
2555
2556  fprintf(treeFile, "\t\"version\": 2,\n");
2557  fprintf(treeFile, "\t\"fields\": [\n");
2558  fprintf(treeFile, "\t\"edge_num\", \"parsimony\"\n");
2559  fprintf(treeFile, "\t]\n");
2560  fprintf(treeFile, "}\n");
2561 
2562  fclose(treeFile);
2563
2564  /* JSON format end */ 
2565           
2566  likelihoodWeightsFile = myfopen(likelihoodWeightsFileName, "wb");
2567 
2568       
2569  for(int i = 0; i < tr->numberOfTipsForInsertion; i++)
2570    {
2571      unsigned int       
2572        lmax;
2573       
2574      int 
2575        validEntries = tr->numberOfBranches;
2576       
2577      for(j =  0; j < tr->numberOfBranches; j++) 
2578        {
2579          inf[j].parsimonyScore = tr->bInf[j].epa->parsimonyScore[i];
2580          inf[j].number = j;
2581        }
2582       
2583      qsort(inf, tr->numberOfBranches, sizeof(infoMP), infoCompare);                 
2584               
2585      j = 0;
2586     
2587      lmax = inf[0].parsimonyScore;
2588     
2589      while(j < validEntries && inf[j].parsimonyScore == lmax)   
2590        {                   
2591          fprintf(likelihoodWeightsFile, "%s I%d %u\n", tr->nameList[tr->inserts[i]], inf[j].number, inf[j].parsimonyScore);
2592          tr->bInf[inf[j].number].epa->countThem[i] = 1;
2593          j++;
2594        }                                 
2595    }
2596     
2597  rax_free(inf);
2598   
2599  fclose(likelihoodWeightsFile); 
2600   
2601 
2602  Tree2StringClassify(tr->tree_string, tr, tr->inserts, FALSE, FALSE, FALSE);
2603  treeFile = fopen(labelledTreeFileName, "wb");
2604  fprintf(treeFile, "%s\n", tr->tree_string);
2605  fclose(treeFile);
2606 
2607
2608  printBothOpen("Equally parsimonious read placements written to file: %s\n\n", likelihoodWeightsFileName);
2609  printBothOpen("Labelled reference tree with branch labels (without query sequences) written to file: %s\n\n", originalLabelledTreeFileName); 
2610  printBothOpen("Labelled reference tree with branch labels in portable pplacer/EPA format (without query sequences) written to file: %s\n\n", jointFormatTreeFileName);
2611  printBothOpen("Labelled reference tree including branch labels and query sequences written to file: %s\n\n", labelledTreeFileName); 
2612 
2613  exit(0);
2614}
2615
2616#ifdef __AVX
2617
2618#ifdef _SSE3_WAS_DEFINED
2619
2620#define __SIM_SSE3
2621
2622#undef _SSE3_WAS_DEFINED
2623
2624#endif
2625
2626#endif
Note: See TracBrowser for help on using the repository browser.