source: branches/nameserver/GDE/RAxML/fastDNAparsimony.c

Last change on this file was 13845, checked in by westram, 10 years ago
  • remove executable flag
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(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(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(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(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(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(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      int 
1756        j = 0;
1757
1758      unsigned char
1759        *nodesInTree = (unsigned char*)rax_calloc((size_t)(tr->mxtips + 1), sizeof(unsigned char));           
1760       
1761      tr->start = findAnyTipFast(tr->start, tr->rdta->numsp);
1762       
1763      tr->bestParsimony = INT_MAX;
1764
1765      evaluateParsimony(tr, tr->start->back, TRUE);
1766               
1767      assert(tr->start);
1768     
1769      checkSeed(adef);
1770
1771      markNodesInTree(tr->start, tr, nodesInTree);
1772      markNodesInTree(tr->start->back, tr, nodesInTree);
1773       
1774      j = tr->ntips + 1;
1775       
1776      if(tr->grouped)
1777        {
1778          for(i = 1; i <= tr->mxtips; i++)     
1779            {
1780              if(tr->constraintVector[i] == -1) 
1781                {
1782                  perm[j++] = i;               
1783                  tr->constraintVector[i] = -9;
1784                }
1785            }
1786        }
1787      else
1788        {
1789          if(tr->constrained)
1790            { 
1791              for(i = 1; i <= tr->mxtips; i++)
1792                tr->constraintVector[i] = 0;
1793             
1794              for(i = 1; i <= tr->mxtips; i++)
1795                {                 
1796                  if(nodesInTree[i] == 0)             
1797                    perm[j++] = i;
1798                  else
1799                    tr->constraintVector[i] = 1;                   
1800                }
1801            }
1802          else
1803            {
1804              for(i = 1; i <= tr->mxtips; i++)     
1805                if(nodesInTree[i] == 0) 
1806                  perm[j++] = i;         
1807            }
1808        }
1809       
1810      for(i = tr->ntips + 1; i <= tr->mxtips; i++) 
1811        {           
1812          int k, j;
1813         
1814          k =  (int)((double)(tr->mxtips + 1 - i) * randum(&adef->parsimonySeed));
1815         
1816          assert(i + k <= tr->mxtips);
1817          j        = perm[i];
1818          perm[i]     = perm[i + k];
1819          perm[i + k] = j;
1820        }   
1821       
1822      f = tr->start;     
1823
1824      rax_free(nodesInTree);
1825    }
1826  else
1827    {
1828      assert(!tr->constrained);
1829
1830      makePermutationFast(perm, tr->mxtips, adef);
1831     
1832      tr->ntips = 0;   
1833     
1834      tr->nextnode = tr->mxtips + 1;       
1835     
1836      buildSimpleTree(tr, perm[1], perm[2], perm[3]);     
1837     
1838      f = tr->start;
1839    }     
1840 
1841  while(tr->ntips < tr->mxtips) 
1842    {   
1843      nodeptr q;
1844     
1845      tr->bestParsimony = INT_MAX;
1846      nextsp = ++(tr->ntips);             
1847      p = tr->nodep[perm[nextsp]];                 
1848      q = tr->nodep[(tr->nextnode)++];
1849      p->back = q;
1850      q->back = p;
1851       
1852      if(tr->grouped && !full)
1853        {
1854          int 
1855            number = p->back->number;           
1856
1857          tr->constraintVector[number] = -9;
1858        }
1859         
1860      stepwiseAddition(tr, q, f->back);                 
1861     
1862      {
1863        nodeptr   
1864          r = tr->insertNode->back;
1865       
1866        int counter = 4;
1867       
1868        hookupDefault(q->next,       tr->insertNode, tr->numBranches);
1869        hookupDefault(q->next->next, r, tr->numBranches);
1870       
1871        computeTraversalInfoParsimony(q, tr->ti, &counter, tr->mxtips, FALSE);             
1872        tr->ti[0] = counter;
1873       
1874        newviewParsimonyIterativeFast(tr);     
1875      }
1876    }   
1877 
1878  //printf("ADD: %d\n", tr->bestParsimony);
1879 
1880  nodeRectifier(tr);
1881 
1882  if(adef->stepwiseAdditionOnly == FALSE)
1883    {
1884      randomMP = tr->bestParsimony;       
1885     
1886      do
1887        {
1888          startMP = randomMP;
1889          nodeRectifier(tr);
1890          for(i = 1; i <= tr->mxtips + tr->mxtips - 2; i++)
1891            {
1892              rearrangeParsimony(tr, tr->nodep[i], 1, 20, FALSE);
1893              if(tr->bestParsimony < randomMP)
1894                {               
1895                  restoreTreeRearrangeParsimony(tr);
1896                  randomMP = tr->bestParsimony;
1897                }
1898            }                             
1899        }
1900      while(randomMP < startMP);
1901    }
1902 
1903  //printf("OPT: %d\n", tr->bestParsimony);
1904
1905 
1906     
1907  rax_free(perm); 
1908  rax_free(tr->parsimonyScore);
1909 
1910  for(model = 0; model < (size_t) tr->NumberOfModels; model++)
1911    rax_free(tr->partitionData[model].parsVect);
1912 
1913  rax_free(tr->ti);
1914} 
1915
1916
1917
1918 
1919static void insertRandom (nodeptr p, nodeptr q, int numBranches)
1920{
1921  nodeptr  r;
1922 
1923  r = q->back;
1924 
1925  hookupDefault(p->next,       q, numBranches);
1926  hookupDefault(p->next->next, r, numBranches); 
1927} 
1928
1929
1930
1931
1932
1933
1934
1935static void buildSimpleTreeRandom (tree *tr, int ip, int iq, int ir)
1936{   
1937  nodeptr  p, s;
1938  int  i;
1939 
1940  i = MIN(ip, iq);
1941  if (ir < i)  i = ir; 
1942  tr->start = tr->nodep[i];
1943  tr->ntips = 3;
1944  p = tr->nodep[ip];
1945  hookupDefault(p, tr->nodep[iq], tr->numBranches);
1946  s = buildNewTip(tr, tr->nodep[ir]);
1947  insertRandom(s, p, tr->numBranches);
1948}
1949
1950int checker(tree *tr, nodeptr p)
1951{
1952  int group = tr->constraintVector[p->number];
1953
1954  if(isTip(p->number, tr->mxtips))
1955    {
1956      group = tr->constraintVector[p->number];
1957      return group;
1958    }
1959  else
1960    {
1961      if(group != -9) 
1962        return group;
1963
1964      group = checker(tr, p->next->back);
1965      if(group != -9) 
1966        return group;
1967
1968      group = checker(tr, p->next->next->back);
1969      if(group != -9) 
1970        return group;
1971
1972      return -9;
1973    }
1974}
1975
1976
1977
1978
1979
1980
1981
1982static int markBranches(nodeptr *branches, nodeptr p, int *counter, int numsp)
1983{
1984  if(isTip(p->number, numsp))
1985    return 0;
1986  else
1987    {
1988      branches[*counter] = p->next;
1989      branches[*counter + 1] = p->next->next;
1990     
1991      *counter = *counter + 2;
1992     
1993      return ((2 + markBranches(branches, p->next->back, counter, numsp) + 
1994               markBranches(branches, p->next->next->back, counter, numsp)));
1995    }
1996}
1997
1998
1999
2000nodeptr findAnyTip(nodeptr p, int numsp)
2001{ 
2002  return  isTip(p->number, numsp) ? p : findAnyTip(p->next->back, numsp);
2003} 
2004
2005
2006int randomInt(int n)
2007{
2008  return rand() %n;
2009}
2010
2011void makePermutation(int *perm, int n, analdef *adef)
2012{   
2013  int  i, j, k;
2014
2015  checkSeed(adef);         
2016
2017  for (i = 1; i <= n; i++)   
2018    perm[i] = i;               
2019
2020  for (i = 1; i <= n; i++) 
2021    {   
2022      k =  (int)((double)(n + 1 - i) * randum(&adef->parsimonySeed));
2023
2024      assert(i + k <= n);
2025     
2026      j        = perm[i];
2027      perm[i]     = perm[i + k];
2028      perm[i + k] = j; 
2029    }
2030}
2031
2032
2033
2034
2035
2036
2037
2038
2039boolean tipHomogeneityChecker(tree *tr, nodeptr p, int grouping)
2040{
2041  if(isTip(p->number, tr->mxtips))
2042    {
2043      if(tr->constraintVector[p->number] != grouping) 
2044        return FALSE;
2045      else 
2046        return TRUE;
2047    }
2048  else
2049    {   
2050      return  (tipHomogeneityChecker(tr, p->next->back, grouping) && tipHomogeneityChecker(tr, p->next->next->back,grouping));     
2051    }
2052}
2053
2054
2055
2056
2057
2058
2059
2060void makeRandomTree(tree *tr, analdef *adef)
2061{ 
2062  nodeptr p, f, randomBranch;   
2063  int nextsp;
2064  int *perm, branchCounter;
2065  nodeptr *branches;
2066 
2067  branches = (nodeptr *)rax_malloc(sizeof(nodeptr) * (2 * tr->mxtips));
2068  perm = (int *)rax_malloc((tr->mxtips + 1) * sizeof(int));                         
2069 
2070  makePermutation(perm, tr->mxtips, adef);             
2071 
2072  tr->ntips = 0;               
2073  tr->nextnode = tr->mxtips + 1;   
2074 
2075  buildSimpleTreeRandom(tr, perm[1], perm[2], perm[3]);
2076 
2077  while (tr->ntips < tr->mxtips) 
2078    {         
2079      tr->bestParsimony = INT_MAX;
2080      nextsp = ++(tr->ntips);             
2081      p = tr->nodep[perm[nextsp]];
2082     
2083      /*printf("ADDING SPECIES %d\n", nextsp);*/
2084     
2085      buildNewTip(tr, p);       
2086     
2087      f = findAnyTip(tr->start, tr->mxtips);
2088      f = f->back;
2089     
2090      branchCounter = 1;
2091      branches[0] = f;
2092      markBranches(branches, f, &branchCounter, tr->mxtips);
2093
2094      assert(branchCounter == ((2 * (tr->ntips - 1)) - 3));
2095     
2096      randomBranch = branches[randomInt(branchCounter)];
2097     
2098      insertRandom(p->back, randomBranch, tr->numBranches);
2099     
2100    }
2101  rax_free(perm);           
2102  rax_free(branches);
2103}
2104
2105
2106
2107void nodeRectifier(tree *tr)
2108{
2109  nodeptr *np = (nodeptr *)rax_malloc(2 * tr->mxtips * sizeof(nodeptr));
2110  int i;
2111  int count = 0;
2112 
2113  tr->start       = tr->nodep[1];
2114  tr->rooted      = FALSE;
2115
2116  /* TODO why is tr->rooted set to FALSE here ?*/
2117 
2118  for(i = tr->mxtips + 1; i <= (tr->mxtips + tr->mxtips - 1); i++)
2119    np[i] = tr->nodep[i];           
2120 
2121  reorderNodes(tr, np, tr->start->back, &count); 
2122
2123 
2124  rax_free(np);
2125}
2126
2127
2128void makeParsimonyTree(tree *tr, analdef *adef)
2129{ 
2130  makeParsimonyTreeFast(tr, adef, TRUE);       
2131}
2132
2133void makeParsimonyTreeIncomplete(tree *tr, analdef *adef)
2134{   
2135  makeParsimonyTreeFast(tr, adef, FALSE);       
2136}
2137 
2138
2139static void setupBranchMetaInfo(tree *tr, nodeptr p, int nTips, branchInfo *bInf)
2140{
2141  int 
2142    countBranches = tr->branchCounter;
2143
2144  if(isTip(p->number, tr->mxtips))   
2145    {     
2146      p->bInf       = &bInf[countBranches];
2147      p->back->bInf = &bInf[countBranches];                           
2148
2149      bInf[countBranches].oP = p;
2150      bInf[countBranches].oQ = p->back;
2151     
2152      bInf[countBranches].epa->leftNodeNumber = p->number;
2153      bInf[countBranches].epa->rightNodeNumber = p->back->number;
2154         
2155      bInf[countBranches].epa->branchNumber = countBranches;                     
2156      bInf[countBranches].epa->originalBranchLength = p->z[0];
2157
2158      tr->branchCounter =  tr->branchCounter + 1;
2159      return;
2160    }
2161  else
2162    {
2163      nodeptr q;
2164      assert(p == p->next->next->next);
2165
2166      p->bInf       = &bInf[countBranches];
2167      p->back->bInf = &bInf[countBranches];
2168
2169      bInf[countBranches].oP = p;
2170      bInf[countBranches].oQ = p->back;
2171
2172      bInf[countBranches].epa->leftNodeNumber = p->number;
2173      bInf[countBranches].epa->rightNodeNumber = p->back->number;
2174
2175               
2176      bInf[countBranches].epa->branchNumber = countBranches;
2177      bInf[countBranches].epa->originalBranchLength = p->z[0];
2178
2179      tr->branchCounter =  tr->branchCounter + 1;     
2180
2181      q = p->next;
2182
2183      while(q != p)
2184        {
2185          setupBranchMetaInfo(tr, q->back, nTips, bInf);       
2186          q = q->next;
2187        }
2188     
2189      return;
2190    }
2191}
2192 
2193
2194
2195static void setupJointFormat(tree *tr, nodeptr p, int ntips, branchInfo *bInf, int *count)
2196{
2197  if(isTip(p->number, tr->mxtips))   
2198    {     
2199      p->bInf->epa->jointLabel = *count;
2200      *count = *count + 1;
2201           
2202      return;
2203    }
2204  else
2205    {                           
2206      setupJointFormat(tr, p->next->back, ntips, bInf, count);           
2207      setupJointFormat(tr, p->next->next->back, ntips, bInf, count);     
2208     
2209      p->bInf->epa->jointLabel = *count;
2210      *count = *count + 1; 
2211     
2212      return;
2213    }
2214}
2215 
2216
2217
2218
2219
2220
2221static void setupBranchInfo(tree *tr, nodeptr q)
2222{
2223  nodeptr
2224    originalNode = tr->nodep[tr->mxtips + 1];
2225
2226  int 
2227    count = 0;
2228
2229  tr->branchCounter = 0;
2230
2231  setupBranchMetaInfo(tr, q, tr->ntips, tr->bInf);
2232   
2233  assert(tr->branchCounter == tr->numberOfBranches);
2234
2235  if(tr->wasRooted)
2236    {
2237      assert(tr->leftRootNode->back == tr->rightRootNode);
2238      assert(tr->leftRootNode       == tr->rightRootNode->back);     
2239
2240      if(!isTip(tr->leftRootNode->number, tr->mxtips))
2241        {
2242          setupJointFormat(tr,  tr->leftRootNode->next->back, tr->ntips, tr->bInf, &count);
2243          setupJointFormat(tr,  tr->leftRootNode->next->next->back, tr->ntips, tr->bInf, &count);
2244        }
2245     
2246       tr->leftRootNode->bInf->epa->jointLabel = count;
2247       tr->rootLabel = count;
2248       count = count + 1;
2249
2250       if(!isTip(tr->rightRootNode->number, tr->mxtips))
2251         {
2252          setupJointFormat(tr,  tr->rightRootNode->next->back, tr->ntips, tr->bInf, &count);
2253          setupJointFormat(tr,  tr->rightRootNode->next->next->back, tr->ntips, tr->bInf, &count);
2254        }             
2255    }
2256  else
2257    {
2258      setupJointFormat(tr, originalNode->back, tr->ntips, tr->bInf, &count);
2259      setupJointFormat(tr, originalNode->next->back, tr->ntips, tr->bInf, &count);
2260      setupJointFormat(tr, originalNode->next->next->back, tr->ntips, tr->bInf, &count);     
2261    } 
2262
2263  assert(count == tr->numberOfBranches);
2264}
2265
2266static void testInsertFast(tree *tr, nodeptr r, nodeptr q)
2267{
2268  unsigned int
2269    result;
2270 
2271  nodeptr 
2272    x = q->back;     
2273 
2274  int 
2275    i,
2276    *inserts = tr->inserts;
2277                   
2278  assert(!tr->grouped);                             
2279 
2280  hookupDefault(r->next,       q, tr->numBranches);
2281  hookupDefault(r->next->next, x, tr->numBranches);                             
2282   
2283  newviewParsimony(tr, r);   
2284   
2285  for(i = 0; i < tr->numberOfTipsForInsertion; i++)
2286    {                               
2287      hookupDefault(r, tr->nodep[inserts[i]], tr->numBranches);
2288     
2289      tr->bestParsimony = INT_MAX;
2290
2291      result = evaluateParsimony(tr, r, FALSE);           
2292     
2293      r->back = (nodeptr) NULL;
2294      tr->nodep[inserts[i]]->back = (nodeptr) NULL;
2295     
2296      tr->bInf[q->bInf->epa->branchNumber].epa->parsimonyScore[i] = result;                     
2297    }
2298 
2299  hookupDefault(q, x, tr->numBranches);
2300 
2301  r->next->next->back = r->next->back = (nodeptr) NULL;
2302 
2303}
2304
2305
2306static void traverseTree(tree *tr, nodeptr r, nodeptr q)
2307{       
2308  testInsertFast(tr, r, q);
2309
2310  if(!isTip(q->number, tr->rdta->numsp))
2311    {   
2312      nodeptr
2313        a = q->next;
2314
2315      while(a != q)
2316        {
2317          traverseTree(tr, r, a->back);
2318          a = a->next;
2319        }     
2320    }
2321} 
2322
2323
2324
2325typedef struct
2326  {
2327    unsigned int parsimonyScore; 
2328    int number;
2329  }
2330  infoMP;
2331
2332
2333static int infoCompare(const void *p1, const void *p2)
2334{
2335  infoMP *rc1 = (infoMP *)p1;
2336  infoMP *rc2 = (infoMP *)p2;
2337
2338  unsigned int i = rc1->parsimonyScore;
2339  unsigned int j = rc2->parsimonyScore;
2340
2341  if (i > j)
2342    return (1);
2343  if (i < j)
2344    return (-1);
2345  return (0);
2346}
2347
2348void classifyMP(tree *tr, analdef *adef)
2349{   
2350  int 
2351    *informative = (int *)rax_malloc(sizeof(int) * (size_t)tr->cdta->endsite),
2352    i, 
2353    j, 
2354    *perm;   
2355 
2356  infoMP
2357    *inf;
2358   
2359  nodeptr     
2360    r, 
2361    q;   
2362
2363  char   
2364    jointFormatTreeFileName[1024], 
2365    originalLabelledTreeFileName[1024],
2366    labelledTreeFileName[1024],
2367    likelihoodWeightsFileName[1024];
2368             
2369  FILE   
2370    *likelihoodWeightsFile,
2371    *treeFile;
2372 
2373  unsigned int 
2374    score;       
2375
2376  assert(adef->restart);
2377
2378  determineUninformativeSites(tr, informative);     
2379
2380  compressDNA(tr, informative, TRUE);
2381
2382  rax_free(informative); 
2383
2384  tr->ti = (int*)rax_malloc(sizeof(int) * 4 * (size_t)tr->mxtips);   
2385 
2386  tr->numberOfBranches = 2 * tr->ntips - 3;
2387
2388  printBothOpen("\nRAxML Evolutionary Placement Algorithm using parsimony\n"); 
2389
2390  tr->bestParsimony = INT_MAX;
2391
2392  score = evaluateParsimony(tr, tr->start->back, TRUE);
2393
2394  printBothOpen("\nparsimony score of reference tree: %u\n\n", score);
2395
2396  perm        = (int *)rax_calloc(((size_t)tr->mxtips) + 1, sizeof(int));
2397  tr->inserts = (int *)rax_calloc((size_t)tr->mxtips,   sizeof(int));
2398
2399  markTips(tr->start,       perm, tr->mxtips);
2400  markTips(tr->start->back, perm ,tr->mxtips);
2401
2402  tr->numberOfTipsForInsertion = 0;
2403
2404  for(i = 1; i <= tr->mxtips; i++)
2405    {
2406      if(perm[i] == 0)
2407        {
2408          tr->inserts[tr->numberOfTipsForInsertion] = i;
2409          tr->numberOfTipsForInsertion = tr->numberOfTipsForInsertion + 1;
2410        }
2411    }   
2412
2413  rax_free(perm);
2414 
2415  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); 
2416
2417  assert(tr->numberOfTipsForInsertion == (tr->mxtips - tr->ntips));     
2418
2419  tr->bInf              = (branchInfo*)rax_malloc(tr->numberOfBranches * sizeof(branchInfo)); 
2420
2421  for(i = 0; i < tr->numberOfBranches; i++)
2422    {     
2423      tr->bInf[i].epa = (epaBranchData*)rax_malloc(sizeof(epaBranchData));                   
2424
2425      tr->bInf[i].epa->parsimonyScore = (unsigned int*)rax_malloc(tr->numberOfTipsForInsertion * sizeof(unsigned int));
2426
2427      for(j = 0; j < tr->numberOfTipsForInsertion; j++)
2428        tr->bInf[i].epa->parsimonyScore[j] = INT_MAX;
2429               
2430      tr->bInf[i].epa->branchNumber = i;
2431      tr->bInf[i].epa->countThem = (int*)rax_calloc(tr->numberOfTipsForInsertion, sizeof(int));
2432     
2433      sprintf(tr->bInf[i].epa->branchLabel, "I%d", i);     
2434    } 
2435
2436  r = tr->nodep[(tr->nextnode)++]; 
2437   
2438
2439  q = findAnyTip(tr->start, tr->rdta->numsp);
2440
2441  assert(isTip(q->number, tr->rdta->numsp));
2442  assert(!isTip(q->back->number, tr->rdta->numsp));
2443         
2444  q = q->back; 
2445 
2446  setupBranchInfo(tr, q);   
2447 
2448  traverseTree(tr, r, q);
2449                   
2450  printBothOpen("Overall Classification time: %f\n\n", gettime() - masterTime); 
2451 
2452 
2453  strcpy(jointFormatTreeFileName,      workdir); 
2454  strcpy(originalLabelledTreeFileName, workdir); 
2455  strcpy(labelledTreeFileName,         workdir); 
2456  strcpy(likelihoodWeightsFileName,    workdir);
2457
2458  strcat(jointFormatTreeFileName,      "RAxML_portableTree."); 
2459  strcat(originalLabelledTreeFileName, "RAxML_originalLabelledTree.");
2460  strcat(labelledTreeFileName,         "RAxML_labelledTree.");
2461  strcat(likelihoodWeightsFileName,    "RAxML_equallyParsimoniousPlacements.");
2462 
2463  strcat(jointFormatTreeFileName,      run_id); 
2464  strcat(originalLabelledTreeFileName, run_id);
2465  strcat(labelledTreeFileName,         run_id);
2466  strcat(likelihoodWeightsFileName,    run_id);
2467
2468  strcat(jointFormatTreeFileName,      ".jplace");
2469 
2470  rax_free(tr->tree_string);
2471
2472  tr->treeStringLength *= 16;
2473
2474  tr->tree_string  = (char*)rax_calloc(tr->treeStringLength, sizeof(char)); 
2475 
2476 
2477
2478  treeFile = myfopen(originalLabelledTreeFileName, "wb");
2479  Tree2StringClassify(tr->tree_string, tr, tr->inserts, TRUE, FALSE, FALSE);
2480  fprintf(treeFile, "%s\n", tr->tree_string);   
2481  fclose(treeFile); 
2482
2483  treeFile = myfopen(jointFormatTreeFileName, "wb");
2484  Tree2StringClassify(tr->tree_string, tr, tr->inserts, TRUE, TRUE, FALSE);
2485 
2486  fprintf(treeFile, "{\n");
2487  fprintf(treeFile, "\t\"tree\": \"%s\", \n", tr->tree_string);
2488  fprintf(treeFile, "\t\"placements\": [\n");
2489     
2490
2491  inf = (infoMP*)rax_malloc(sizeof(infoMP) * tr->numberOfBranches); 
2492 
2493  /* joint format */
2494       
2495 
2496
2497  for(i = 0; i < tr->numberOfTipsForInsertion; i++)   
2498    {
2499      unsigned int
2500        lmax;
2501     
2502      int         
2503        validEntries = tr->numberOfBranches;
2504     
2505      for(j =  0; j < tr->numberOfBranches; j++) 
2506        {
2507          inf[j].parsimonyScore = tr->bInf[j].epa->parsimonyScore[i];     
2508          inf[j].number         = tr->bInf[j].epa->jointLabel;
2509        }
2510     
2511      qsort(inf, tr->numberOfBranches, sizeof(infoMP), infoCompare);               
2512     
2513      j = 0;
2514     
2515      lmax = inf[0].parsimonyScore;
2516     
2517      fprintf(treeFile, "\t{\"p\":[");
2518     
2519      while(j < validEntries && inf[j].parsimonyScore == lmax)   
2520        {           
2521          if(j > 0)
2522            {
2523              if(tr->wasRooted && inf[j].number == tr->rootLabel)                 
2524                assert(0);               
2525              else
2526                fprintf(treeFile, ",[%d, %u]", inf[j].number, inf[j].parsimonyScore);
2527            }
2528          else
2529            {
2530              if(tr->wasRooted && inf[j].number == tr->rootLabel)                 
2531                assert(0);               
2532              else
2533                fprintf(treeFile, "[%d, %u]", inf[j].number, inf[j].parsimonyScore);
2534            }     
2535         
2536          j++;
2537        }
2538     
2539      if(i == tr->numberOfTipsForInsertion - 1)
2540        fprintf(treeFile, "], \"n\":[\"%s\"]}\n", tr->nameList[tr->inserts[i]]);
2541      else
2542        fprintf(treeFile, "], \"n\":[\"%s\"]},\n", tr->nameList[tr->inserts[i]]);
2543    } 
2544   
2545  fprintf(treeFile, "\t ],\n");
2546  /*  fprintf(treeFile, "\t\"metadata\": {\"invocation\": \"RAxML EPA parsimony\"},\n");*/
2547  fprintf(treeFile, "\t\"metadata\": {\"invocation\": ");
2548
2549  fprintf(treeFile, "\"");
2550 
2551  {
2552    int i;
2553   
2554    for(i = 0; i < globalArgc; i++)
2555      fprintf(treeFile,"%s ", globalArgv[i]);
2556  }
2557  fprintf(treeFile, "\", \"raxml_version\": \"%s\"", programVersion);
2558  fprintf(treeFile,"},\n");
2559
2560  fprintf(treeFile, "\t\"version\": 2,\n");
2561  fprintf(treeFile, "\t\"fields\": [\n");
2562  fprintf(treeFile, "\t\"edge_num\", \"parsimony\"\n");
2563  fprintf(treeFile, "\t]\n");
2564  fprintf(treeFile, "}\n");
2565 
2566  fclose(treeFile);
2567
2568  /* JSON format end */ 
2569           
2570  likelihoodWeightsFile = myfopen(likelihoodWeightsFileName, "wb");
2571 
2572       
2573  for(i = 0; i < tr->numberOfTipsForInsertion; i++)   
2574    {
2575      unsigned int       
2576        lmax;
2577       
2578      int 
2579        validEntries = tr->numberOfBranches;
2580       
2581      for(j =  0; j < tr->numberOfBranches; j++) 
2582        {
2583          inf[j].parsimonyScore = tr->bInf[j].epa->parsimonyScore[i];
2584          inf[j].number = j;
2585        }
2586       
2587      qsort(inf, tr->numberOfBranches, sizeof(infoMP), infoCompare);                 
2588               
2589      j = 0;
2590     
2591      lmax = inf[0].parsimonyScore;
2592     
2593      while(j < validEntries && inf[j].parsimonyScore == lmax)   
2594        {                   
2595          fprintf(likelihoodWeightsFile, "%s I%d %u\n", tr->nameList[tr->inserts[i]], inf[j].number, inf[j].parsimonyScore);
2596          tr->bInf[inf[j].number].epa->countThem[i] = 1;
2597          j++;
2598        }                                 
2599    }
2600     
2601  rax_free(inf);
2602   
2603  fclose(likelihoodWeightsFile); 
2604   
2605 
2606  Tree2StringClassify(tr->tree_string, tr, tr->inserts, FALSE, FALSE, FALSE);
2607  treeFile = fopen(labelledTreeFileName, "wb");
2608  fprintf(treeFile, "%s\n", tr->tree_string);
2609  fclose(treeFile);
2610 
2611
2612  printBothOpen("Equally parsimonious read placements written to file: %s\n\n", likelihoodWeightsFileName);
2613  printBothOpen("Labelled reference tree with branch labels (without query sequences) written to file: %s\n\n", originalLabelledTreeFileName); 
2614  printBothOpen("Labelled reference tree with branch labels in portable pplacer/EPA format (without query sequences) written to file: %s\n\n", jointFormatTreeFileName);
2615  printBothOpen("Labelled reference tree including branch labels and query sequences written to file: %s\n\n", labelledTreeFileName); 
2616 
2617  exit(0);
2618}
2619
2620#ifdef __AVX
2621
2622#ifdef _SSE3_WAS_DEFINED
2623
2624#define __SIM_SSE3
2625
2626#undef _SSE3_WAS_DEFINED
2627
2628#endif
2629
2630#endif
Note: See TracBrowser for help on using the repository browser.