1
/* RAxML-VI-HPC (version 2.2) a program for sequential and parallel estimation of phylogenetic trees
2
* Copyright August 2006 by Alexandros Stamatakis
4
* Partially derived from
5
* fastDNAml, a program for estimation of phylogenetic trees from sequences by Gary J. Olsen
9
* Programs of the PHYLIP package by Joe Felsenstein.
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)
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
22
* For any other enquiries send an Email to Alexandros Stamatakis
23
* Alexandros.Stamatakis@epfl.ch
25
* When publishing work that is based on the results from RAxML-VI-HPC please cite:
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
34
#include <sys/times.h>
35
#include <sys/types.h>
53
#define _SSE3_WAS_DEFINED
64
#include <xmmintrin.h>
65
#include <pmmintrin.h>
71
#include <xmmintrin.h>
72
#include <immintrin.h>
79
extern const unsigned int mask32[32];
80
/* vector-specific stuff */
82
extern char **globalArgv;
83
extern int globalArgc;
87
#define INTS_PER_VECTOR 4
88
#define INT_TYPE __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
102
#define INTS_PER_VECTOR 8
103
#define INT_TYPE __m256d
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
115
extern double masterTime;
116
extern char workdir[1024];
117
extern char run_id[128];
119
/********************************DNA FUNCTIONS *****************************************************************/
123
static void checkSeed(analdef *adef)
125
static boolean seedChecked = FALSE;
129
/*printf("Checking seed\n");*/
131
if(adef->parsimonySeed <= 0)
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");
137
assert(adef->parsimonySeed > 0);
142
static void getxnodeLocal (nodeptr p)
146
if((s = p->next)->x || (s = s->next)->x)
153
static void computeTraversalInfoParsimony(nodeptr p, int *ti, int *counter, int maxTips, boolean full)
157
r = p->next->next->back;
164
if(q->number > maxTips)
165
computeTraversalInfoParsimony(q, ti, counter, maxTips, full);
167
if(r->number > maxTips)
168
computeTraversalInfoParsimony(r, ti, counter, maxTips, full);
172
if(q->number > maxTips && !q->x)
173
computeTraversalInfoParsimony(q, ti, counter, maxTips, full);
175
if(r->number > maxTips && !r->x)
176
computeTraversalInfoParsimony(r, ti, counter, maxTips, full);
180
ti[*counter] = p->number;
181
ti[*counter + 1] = q->number;
182
ti[*counter + 2] = r->number;
183
*counter = *counter + 4;
188
#if (defined(__SIM_SSE3) || defined(__AVX))
190
static inline unsigned int populationCount(INT_TYPE v_N)
193
res[INTS_PER_VECTOR] __attribute__ ((aligned (BYTE_ALIGNMENT)));
199
VECTOR_STORE((CAST)res, v_N);
201
for(i = 0; i < INTS_PER_VECTOR; i++)
202
a += BIT_COUNT(res[i]);
209
static inline unsigned int populationCount(unsigned int n)
216
#if (defined(__SIM_SSE3) || defined(__AVX))
218
void newviewParsimonyIterativeFast(tree *tr)
221
allOne = SET_ALL_BITS_ONE;
229
for(index = 4; index < count; index += 4)
235
pNumber = (size_t)ti[index],
236
qNumber = (size_t)ti[index + 1],
237
rNumber = (size_t)ti[index + 2];
239
for(model = 0; model < tr->NumberOfModels; model++)
243
states = tr->partitionData[model].states,
244
width = tr->partitionData[model].parsimonyLength;
258
for(k = 0; k < 2; k++)
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]);
265
for(i = 0; i < width; i += INTS_PER_VECTOR)
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);
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);
282
v_N = VECTOR_BIT_OR(l_A, l_C);
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)));
287
v_N = VECTOR_AND_NOT(v_N, allOne);
289
totalScore += populationCount(v_N);
300
for(k = 0; k < 4; k++)
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]);
307
for(i = 0; i < width; i += INTS_PER_VECTOR)
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);
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);
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);
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);
334
v_N = VECTOR_BIT_OR(VECTOR_BIT_OR(l_A, l_C), VECTOR_BIT_OR(l_G, l_T));
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)));
341
v_N = VECTOR_AND_NOT(v_N, allOne);
343
totalScore += populationCount(v_N);
354
for(k = 0; k < 20; k++)
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]);
361
for(i = 0; i < width; i += INTS_PER_VECTOR)
367
v_N = SET_ALL_BITS_ZERO,
371
for(j = 0; j < 20; j++)
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);
378
v_N = VECTOR_BIT_OR(v_N, l_A[j]);
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])));
384
v_N = VECTOR_AND_NOT(v_N, allOne);
386
totalScore += populationCount(v_N);
397
assert(states <= 32);
399
for(k = 0; k < states; k++)
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]);
406
for(i = 0; i < width; i += INTS_PER_VECTOR)
412
v_N = SET_ALL_BITS_ZERO,
416
for(j = 0; j < states; j++)
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);
423
v_N = VECTOR_BIT_OR(v_N, l_A[j]);
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])));
429
v_N = VECTOR_AND_NOT(v_N, allOne);
431
totalScore += populationCount(v_N);
437
tr->parsimonyScore[pNumber] = totalScore + tr->parsimonyScore[rNumber] + tr->parsimonyScore[qNumber];
443
unsigned int evaluateParsimonyIterativeFast(tree *tr)
446
allOne = SET_ALL_BITS_ONE;
449
pNumber = (size_t)tr->ti[1],
450
qNumber = (size_t)tr->ti[2];
456
bestScore = tr->bestParsimony,
460
newviewParsimonyIterativeFast(tr);
462
sum = tr->parsimonyScore[pNumber] + tr->parsimonyScore[qNumber];
464
for(model = 0; model < tr->NumberOfModels; model++)
468
states = tr->partitionData[model].states,
469
width = tr->partitionData[model].parsimonyLength,
480
for(k = 0; k < 2; k++)
482
left[k] = &(tr->partitionData[model].parsVect[(width * 2 * qNumber) + width * k]);
483
right[k] = &(tr->partitionData[model].parsVect[(width * 2 * pNumber) + width * k]);
486
for(i = 0; i < width; i += INTS_PER_VECTOR)
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);
493
v_N = VECTOR_AND_NOT(v_N, allOne);
495
sum += populationCount(v_N);
508
for(k = 0; k < 4; k++)
510
left[k] = &(tr->partitionData[model].parsVect[(width * 4 * qNumber) + width * k]);
511
right[k] = &(tr->partitionData[model].parsVect[(width * 4 * pNumber) + width * k]);
514
for(i = 0; i < width; i += INTS_PER_VECTOR)
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));
523
v_N = VECTOR_AND_NOT(v_N, allOne);
525
sum += populationCount(v_N);
538
for(k = 0; k < 20; k++)
540
left[k] = &(tr->partitionData[model].parsVect[(width * 20 * qNumber) + width * k]);
541
right[k] = &(tr->partitionData[model].parsVect[(width * 20 * pNumber) + width * k]);
544
for(i = 0; i < width; i += INTS_PER_VECTOR)
551
v_N = SET_ALL_BITS_ZERO;
553
for(j = 0; j < 20; j++)
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);
559
v_N = VECTOR_AND_NOT(v_N, allOne);
561
sum += populationCount(v_N);
574
assert(states <= 32);
576
for(k = 0; k < states; k++)
578
left[k] = &(tr->partitionData[model].parsVect[(width * states * qNumber) + width * k]);
579
right[k] = &(tr->partitionData[model].parsVect[(width * states * pNumber) + width * k]);
582
for(i = 0; i < width; i += INTS_PER_VECTOR)
589
v_N = SET_ALL_BITS_ZERO;
591
for(j = 0; j < states; j++)
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);
597
v_N = VECTOR_AND_NOT(v_N, allOne);
599
sum += populationCount(v_N);
614
void newviewParsimonyIterativeFast(tree *tr)
622
for(index = 4; index < count; index += 4)
628
pNumber = (size_t)ti[index],
629
qNumber = (size_t)ti[index + 1],
630
rNumber = (size_t)ti[index + 2];
632
for(model = 0; model < tr->NumberOfModels; model++)
636
states = tr->partitionData[model].states,
637
width = tr->partitionData[model].parsimonyLength;
658
for(k = 0; k < 2; k++)
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]);
665
for(i = 0; i < width; i++)
667
t_A = left[0][i] & right[0][i];
668
t_C = left[1][i] & right[1][i];
670
o_A = left[0][i] | right[0][i];
671
o_C = left[1][i] | right[1][i];
675
this[0][i] = t_A | (t_N & o_A);
676
this[1][i] = t_C | (t_N & o_C);
678
totalScore += populationCount(t_N);
689
for(k = 0; k < 4; k++)
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]);
707
for(i = 0; i < width; i++)
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];
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];
719
t_N = ~(t_A | t_C | t_G | t_T);
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);
726
totalScore += populationCount(t_N);
742
for(k = 0; k < 20; k++)
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]);
749
for(i = 0; i < width; i++)
755
for(k = 0; k < 20; k++)
757
t_A[k] = left[k][i] & right[k][i];
758
o_A[k] = left[k][i] | right[k][i];
764
for(k = 0; k < 20; k++)
765
this[k][i] = t_A[k] | (t_N & o_A[k]);
767
totalScore += populationCount(t_N);
783
assert(states <= 32);
785
for(k = 0; k < states; k++)
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]);
792
for(i = 0; i < width; i++)
796
for(k = 0; k < states; k++)
798
t_A[k] = left[k][i] & right[k][i];
799
o_A[k] = left[k][i] | right[k][i];
805
for(k = 0; k < states; k++)
806
this[k][i] = t_A[k] | (t_N & o_A[k]);
808
totalScore += populationCount(t_N);
814
tr->parsimonyScore[pNumber] = totalScore + tr->parsimonyScore[rNumber] + tr->parsimonyScore[qNumber];
820
unsigned int evaluateParsimonyIterativeFast(tree *tr)
823
pNumber = (size_t)tr->ti[1],
824
qNumber = (size_t)tr->ti[2];
830
bestScore = tr->bestParsimony,
834
newviewParsimonyIterativeFast(tr);
836
sum = tr->parsimonyScore[pNumber] + tr->parsimonyScore[qNumber];
838
for(model = 0; model < tr->NumberOfModels; model++)
842
states = tr->partitionData[model].states,
843
width = tr->partitionData[model].parsimonyLength,
857
for(k = 0; k < 2; k++)
859
left[k] = &(tr->partitionData[model].parsVect[(width * 2 * qNumber) + width * k]);
860
right[k] = &(tr->partitionData[model].parsVect[(width * 2 * pNumber) + width * k]);
863
for(i = 0; i < width; i++)
865
t_A = left[0][i] & right[0][i];
866
t_C = left[1][i] & right[1][i];
870
sum += populationCount(t_N);
888
for(k = 0; k < 4; k++)
890
left[k] = &(tr->partitionData[model].parsVect[(width * 4 * qNumber) + width * k]);
891
right[k] = &(tr->partitionData[model].parsVect[(width * 4 * pNumber) + width * k]);
894
for(i = 0; i < width; i++)
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];
901
t_N = ~(t_A | t_C | t_G | t_T);
903
sum += populationCount(t_N);
918
for(k = 0; k < 20; k++)
920
left[k] = &(tr->partitionData[model].parsVect[(width * 20 * qNumber) + width * k]);
921
right[k] = &(tr->partitionData[model].parsVect[(width * 20 * pNumber) + width * k]);
924
for(i = 0; i < width; i++)
928
for(k = 0; k < 20; k++)
930
t_A = left[k][i] & right[k][i];
936
sum += populationCount(t_N);
951
assert(states <= 32);
953
for(k = 0; k < states; k++)
955
left[k] = &(tr->partitionData[model].parsVect[(width * states * qNumber) + width * k]);
956
right[k] = &(tr->partitionData[model].parsVect[(width * states * pNumber) + width * k]);
959
for(i = 0; i < width; i++)
963
for(k = 0; k < states; k++)
965
t_A = left[k][i] & right[k][i];
971
sum += populationCount(t_N);
990
static unsigned int evaluateParsimony(tree *tr, nodeptr p, boolean full)
992
volatile unsigned int result;
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);
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);
1018
result = evaluateParsimonyIterativeFast(tr);
1024
static void newviewParsimony(tree *tr, nodeptr p)
1026
if(p->number <= tr->mxtips)
1033
computeTraversalInfoParsimony(p, tr->ti, &counter, tr->mxtips, FALSE);
1034
tr->ti[0] = counter;
1036
newviewParsimonyIterativeFast(tr);
1044
/****************************************************************************************************************************************/
1046
static void insertParsimony (tree *tr, nodeptr p, nodeptr q)
1052
hookupDefault(p->next, q, tr->numBranches);
1053
hookupDefault(p->next->next, r, tr->numBranches);
1055
newviewParsimony(tr, p);
1059
static nodeptr buildNewTip (tree *tr, nodeptr p)
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);
1073
static nodeptr buildNewTip (tree *tr, nodeptr p)
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;
1085
static void buildSimpleTree (tree *tr, int ip, int iq, int ir)
1092
tr->start = tr->nodep[i];
1095
hookupDefault(p, tr->nodep[iq], tr->numBranches);
1096
s = buildNewTip(tr, tr->nodep[ir]);
1097
insertParsimony(tr, s, p);
1101
static void testInsertParsimony (tree *tr, nodeptr p, nodeptr q)
1115
rNumber = tr->constraintVector[r->number],
1116
qNumber = tr->constraintVector[q->number],
1117
pNumber = tr->constraintVector[p->number];
1122
pNumber = checker(tr, p->back);
1128
qNumber = checker(tr, q);
1131
rNumber = checker(tr, r);
1133
if(pNumber == rNumber || pNumber == qNumber)
1140
insertParsimony(tr, p, q);
1142
mp = evaluateParsimony(tr, p->next->next, FALSE);
1144
if(mp < tr->bestParsimony)
1146
tr->bestParsimony = mp;
1151
hookupDefault(q, r, tr->numBranches);
1152
p->next->next->back = p->next->back = (nodeptr) NULL;
1159
static void restoreTreeParsimony(tree *tr, nodeptr p, nodeptr q)
1166
hookupDefault(p->next, q, tr->numBranches);
1167
hookupDefault(p->next->next, r, tr->numBranches);
1169
computeTraversalInfoParsimony(p, tr->ti, &counter, tr->mxtips, FALSE);
1170
tr->ti[0] = counter;
1172
newviewParsimonyIterativeFast(tr);
1176
static void addTraverseParsimony (tree *tr, nodeptr p, nodeptr q, int mintrav, int maxtrav, boolean doAll)
1178
if (doAll || (--mintrav <= 0))
1179
testInsertParsimony(tr, p, q);
1181
if (((q->number > tr->mxtips)) && ((--maxtrav > 0) || doAll))
1183
addTraverseParsimony(tr, p, q->next->back, mintrav, maxtrav, doAll);
1184
addTraverseParsimony(tr, p, q->next->next->back, mintrav, maxtrav, doAll);
1189
static nodeptr findAnyTipFast(nodeptr p, int numsp)
1191
return (p->number <= numsp)? p : findAnyTipFast(p->next->back, numsp);
1195
static void makePermutationFast(int *perm, int n, analdef *adef)
1204
for (i = 1; i <= n; i++)
1207
for (i = 1; i <= n; i++)
1209
double d = randum(&adef->parsimonySeed);
1211
k = (int)((double)(n + 1 - i) * d);
1215
perm[i] = perm[i + k];
1220
static nodeptr removeNodeParsimony (nodeptr p, tree *tr)
1225
r = p->next->next->back;
1227
hookupDefault(q, r, tr->numBranches);
1229
p->next->next->back = p->next->back = (node *) NULL;
1234
static int rearrangeParsimony(tree *tr, nodeptr p, int mintrav, int maxtrav, boolean doAll)
1250
if (maxtrav > tr->ntips - 3)
1251
maxtrav = tr->ntips - 3;
1253
assert(mintrav == 1);
1255
if(maxtrav < mintrav)
1262
if(! tipHomogeneityChecker(tr, p->back, 0))
1265
if(! tipHomogeneityChecker(tr, q->back, 0))
1268
if(doQ == FALSE && doP == FALSE)
1272
if((p->number > tr->mxtips) && doP)
1275
p2 = p->next->next->back;
1277
if ((p1->number > tr->mxtips) || (p2->number > tr->mxtips))
1279
removeNodeParsimony(p, tr);
1281
if ((p1->number > tr->mxtips))
1283
addTraverseParsimony(tr, p, p1->next->back, mintrav, maxtrav, doAll);
1284
addTraverseParsimony(tr, p, p1->next->next->back, mintrav, maxtrav, doAll);
1287
if ((p2->number > tr->mxtips))
1289
addTraverseParsimony(tr, p, p2->next->back, mintrav, maxtrav, doAll);
1290
addTraverseParsimony(tr, p, p2->next->next->back, mintrav, maxtrav, doAll);
1294
hookupDefault(p->next, p1, tr->numBranches);
1295
hookupDefault(p->next->next, p2, tr->numBranches);
1297
newviewParsimony(tr, p);
1301
if ((q->number > tr->mxtips) && (maxtrav > 0) && doQ)
1304
q2 = q->next->next->back;
1308
(q1->number > tr->mxtips) &&
1309
((q1->next->back->number > tr->mxtips) || (q1->next->next->back->number > tr->mxtips))
1313
(q2->number > tr->mxtips) &&
1314
((q2->next->back->number > tr->mxtips) || (q2->next->next->back->number > tr->mxtips))
1319
removeNodeParsimony(q, tr);
1321
mintrav2 = mintrav > 2 ? mintrav : 2;
1323
if ((q1->number > tr->mxtips))
1325
addTraverseParsimony(tr, q, q1->next->back, mintrav2 , maxtrav, doAll);
1326
addTraverseParsimony(tr, q, q1->next->next->back, mintrav2 , maxtrav, doAll);
1329
if ((q2->number > tr->mxtips))
1331
addTraverseParsimony(tr, q, q2->next->back, mintrav2 , maxtrav, doAll);
1332
addTraverseParsimony(tr, q, q2->next->next->back, mintrav2 , maxtrav, doAll);
1335
hookupDefault(q->next, q1, tr->numBranches);
1336
hookupDefault(q->next->next, q2, tr->numBranches);
1338
newviewParsimony(tr, q);
1346
static void restoreTreeRearrangeParsimony(tree *tr)
1348
removeNodeParsimony(tr->removeNode, tr);
1349
restoreTreeParsimony(tr, tr->removeNode, tr->insertNode);
1353
static boolean isInformative2(tree *tr, int site)
1356
informativeCounter = 0,
1365
for(j = 0; j < 256; j++)
1368
for(j = 1; j <= tr->mxtips; j++)
1370
nucleotide = tr->yVector[j][site];
1371
check[nucleotide] = check[nucleotide] + 1;
1377
informativeCounter++;
1378
target = target | 1;
1382
informativeCounter++;
1383
target = target | 2;
1387
informativeCounter++;
1388
target = target | 4;
1392
informativeCounter++;
1393
target = target | 8;
1396
if(informativeCounter >= 2)
1400
for(j = 0; j < undetermined; j++)
1402
if(j == 3 || j == 5 || j == 6 || j == 7 || j == 9 || j == 10 || j == 11 ||
1403
j == 12 || j == 13 || j == 14)
1418
static boolean isInformative(tree *tr, int dataType, int site)
1421
informativeCounter = 0,
1424
undetermined = getUndetermined(dataType);
1427
*bitVector = getBitVector(dataType);
1433
for(j = 0; j < 256; j++)
1436
for(j = 1; j <= tr->mxtips; j++)
1438
nucleotide = tr->yVector[j][site];
1439
check[nucleotide] = check[nucleotide] + 1;
1440
assert(bitVector[nucleotide] > 0);
1443
for(j = 0; j < undetermined; j++)
1446
informativeCounter++;
1449
if(informativeCounter <= 1)
1453
for(j = 0; j < undetermined; j++)
1464
static void determineUninformativeSites(tree *tr, int *informative)
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.
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.
1483
for(i = 0; i < tr->cdta->endsite; i++)
1485
if(isInformative(tr, tr->dataVector[i], i))
1495
/* printf("Uninformative Patterns: %d\n", number); */
1499
static void reorderNodes(tree *tr, nodeptr *np, nodeptr p, int *count)
1503
if((p->number <= tr->mxtips))
1507
for(i = tr->mxtips + 1; (i <= (tr->mxtips + tr->mxtips - 1)) && (found == 0); i++)
1509
if (p == np[i] || p == np[i]->next || p == np[i]->next->next)
1512
tr->nodep[*count + tr->mxtips + 1] = np[i];
1515
if(p == np[i]->next)
1516
tr->nodep[*count + tr->mxtips + 1] = np[i]->next;
1518
tr->nodep[*count + tr->mxtips + 1] = np[i]->next->next;
1522
*count = *count + 1;
1528
reorderNodes(tr, np, p->next->back, count);
1529
reorderNodes(tr, np, p->next->next->back, count);
1538
static void compressDNA(tree *tr, int *informative, boolean saveMemory)
1546
totalNodes = (size_t)tr->innerNodes + 1 + (size_t)tr->mxtips;
1548
totalNodes = 2 * (size_t)tr->mxtips;
1552
for(model = 0; model < (size_t) tr->NumberOfModels; model++)
1556
states = (size_t)tr->partitionData[model].states,
1558
compressedEntriesPadded,
1560
lower = tr->partitionData[model].lower,
1561
upper = tr->partitionData[model].upper;
1564
**compressedTips = (parsimonyNumber **)rax_malloc(states * sizeof(parsimonyNumber*)),
1565
*compressedValues = (parsimonyNumber *)rax_malloc(states * sizeof(parsimonyNumber));
1567
for(i = lower; i < upper; i++)
1569
entries += (size_t)tr->cdta->aliaswgt[i];
1571
compressedEntries = entries / PCF;
1573
if(entries % PCF != 0)
1574
compressedEntries++;
1576
#if (defined(__SIM_SSE3) || defined(__AVX))
1577
if(compressedEntries % INTS_PER_VECTOR != 0)
1578
compressedEntriesPadded = compressedEntries + (INTS_PER_VECTOR - (compressedEntries % INTS_PER_VECTOR));
1580
compressedEntriesPadded = compressedEntries;
1582
compressedEntriesPadded = compressedEntries;
1586
tr->partitionData[model].parsVect = (parsimonyNumber *)rax_malloc((size_t)compressedEntriesPadded * states * totalNodes * sizeof(parsimonyNumber));
1588
for(i = 0; i < compressedEntriesPadded * states * totalNodes; i++)
1589
tr->partitionData[model].parsVect[i] = 0;
1591
for(i = 0; i < (size_t)tr->mxtips; i++)
1595
compressedIndex = 0,
1596
compressedCounter = 0,
1599
for(k = 0; k < states; k++)
1601
compressedTips[k] = &(tr->partitionData[model].parsVect[(compressedEntriesPadded * states * (i + 1)) + (compressedEntriesPadded * k)]);
1602
compressedValues[k] = 0;
1605
for(index = lower; index < (size_t)upper; index++)
1607
if(informative[index])
1610
*bitValue = getBitVector(tr->partitionData[model].dataType);
1613
value = bitValue[tr->yVector[i + 1][index]];
1615
for(w = 0; w < (size_t)tr->cdta->aliaswgt[index]; w++)
1617
for(k = 0; k < states; k++)
1619
if(value & mask32[k])
1620
compressedValues[k] |= mask32[compressedCounter];
1623
compressedCounter++;
1625
if(compressedCounter == PCF)
1627
for(k = 0; k < states; k++)
1629
compressedTips[k][compressedIndex] = compressedValues[k];
1630
compressedValues[k] = 0;
1633
compressedCounter = 0;
1640
for(;compressedIndex < compressedEntriesPadded; compressedIndex++)
1642
for(;compressedCounter < PCF; compressedCounter++)
1643
for(k = 0; k < states; k++)
1644
compressedValues[k] |= mask32[compressedCounter];
1646
for(k = 0; k < states; k++)
1648
compressedTips[k][compressedIndex] = compressedValues[k];
1649
compressedValues[k] = 0;
1652
compressedCounter = 0;
1656
tr->partitionData[model].parsimonyLength = compressedEntriesPadded;
1658
rax_free(compressedTips);
1659
rax_free(compressedValues);
1662
tr->parsimonyScore = (unsigned int*)rax_malloc(sizeof(unsigned int) * totalNodes);
1664
for(i = 0; i < totalNodes; i++)
1665
tr->parsimonyScore[i] = 0;
1670
static void stepwiseAddition(tree *tr, nodeptr p, nodeptr q)
1684
p->next->next->back = r;
1685
r->back = p->next->next;
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;
1692
mp = evaluateParsimonyIterativeFast(tr);
1694
if(mp < tr->bestParsimony)
1696
tr->bestParsimony = mp;
1703
if(q->number > tr->mxtips && tr->parsimonyScore[q->number] > 0)
1705
stepwiseAddition(tr, p, q->next->back);
1706
stepwiseAddition(tr, p, q->next->next->back);
1710
static void markNodesInTree(nodeptr p, tree *tr, unsigned char *nodesInTree)
1712
if(isTip(p->number, tr->mxtips))
1713
nodesInTree[p->number] = 1;
1716
markNodesInTree(p->next->back, tr, nodesInTree);
1717
markNodesInTree(p->next->next->back, tr, nodesInTree);
1722
void makeParsimonyTreeFast(tree *tr, analdef *adef, boolean full)
1734
*perm = (int *)rax_malloc((size_t)(tr->mxtips + 1) * sizeof(int)),
1735
*informative = (int *)rax_malloc(sizeof(int) * (size_t)tr->cdta->endsite);
1743
determineUninformativeSites(tr, informative);
1745
compressDNA(tr, informative, FALSE);
1747
rax_free(informative);
1749
tr->ti = (int*)rax_malloc(sizeof(int) * 4 * (size_t)tr->mxtips);
1759
*nodesInTree = (unsigned char*)rax_calloc((size_t)(tr->mxtips + 1), sizeof(unsigned char));
1761
tr->start = findAnyTipFast(tr->start, tr->rdta->numsp);
1763
tr->bestParsimony = INT_MAX;
1765
evaluateParsimony(tr, tr->start->back, TRUE);
1771
markNodesInTree(tr->start, tr, nodesInTree);
1772
markNodesInTree(tr->start->back, tr, nodesInTree);
1778
for(i = 1; i <= tr->mxtips; i++)
1780
if(tr->constraintVector[i] == -1)
1783
tr->constraintVector[i] = -9;
1791
for(i = 1; i <= tr->mxtips; i++)
1792
tr->constraintVector[i] = 0;
1794
for(i = 1; i <= tr->mxtips; i++)
1796
if(nodesInTree[i] == 0)
1799
tr->constraintVector[i] = 1;
1804
for(i = 1; i <= tr->mxtips; i++)
1805
if(nodesInTree[i] == 0)
1810
for(i = tr->ntips + 1; i <= tr->mxtips; i++)
1814
k = (int)((double)(tr->mxtips + 1 - i) * randum(&adef->parsimonySeed));
1816
assert(i + k <= tr->mxtips);
1818
perm[i] = perm[i + k];
1824
rax_free(nodesInTree);
1828
assert(!tr->constrained);
1830
makePermutationFast(perm, tr->mxtips, adef);
1834
tr->nextnode = tr->mxtips + 1;
1836
buildSimpleTree(tr, perm[1], perm[2], perm[3]);
1841
while(tr->ntips < tr->mxtips)
1845
tr->bestParsimony = INT_MAX;
1846
nextsp = ++(tr->ntips);
1847
p = tr->nodep[perm[nextsp]];
1848
q = tr->nodep[(tr->nextnode)++];
1852
if(tr->grouped && !full)
1855
number = p->back->number;
1857
tr->constraintVector[number] = -9;
1860
stepwiseAddition(tr, q, f->back);
1864
r = tr->insertNode->back;
1868
hookupDefault(q->next, tr->insertNode, tr->numBranches);
1869
hookupDefault(q->next->next, r, tr->numBranches);
1871
computeTraversalInfoParsimony(q, tr->ti, &counter, tr->mxtips, FALSE);
1872
tr->ti[0] = counter;
1874
newviewParsimonyIterativeFast(tr);
1878
//printf("ADD: %d\n", tr->bestParsimony);
1882
if(adef->stepwiseAdditionOnly == FALSE)
1884
randomMP = tr->bestParsimony;
1890
for(i = 1; i <= tr->mxtips + tr->mxtips - 2; i++)
1892
rearrangeParsimony(tr, tr->nodep[i], 1, 20, FALSE);
1893
if(tr->bestParsimony < randomMP)
1895
restoreTreeRearrangeParsimony(tr);
1896
randomMP = tr->bestParsimony;
1900
while(randomMP < startMP);
1903
//printf("OPT: %d\n", tr->bestParsimony);
1908
rax_free(tr->parsimonyScore);
1910
for(model = 0; model < (size_t) tr->NumberOfModels; model++)
1911
rax_free(tr->partitionData[model].parsVect);
1919
static void insertRandom (nodeptr p, nodeptr q, int numBranches)
1925
hookupDefault(p->next, q, numBranches);
1926
hookupDefault(p->next->next, r, numBranches);
1935
static void buildSimpleTreeRandom (tree *tr, int ip, int iq, int ir)
1942
tr->start = tr->nodep[i];
1945
hookupDefault(p, tr->nodep[iq], tr->numBranches);
1946
s = buildNewTip(tr, tr->nodep[ir]);
1947
insertRandom(s, p, tr->numBranches);
1950
int checker(tree *tr, nodeptr p)
1952
int group = tr->constraintVector[p->number];
1954
if(isTip(p->number, tr->mxtips))
1956
group = tr->constraintVector[p->number];
1964
group = checker(tr, p->next->back);
1968
group = checker(tr, p->next->next->back);
1982
static int markBranches(nodeptr *branches, nodeptr p, int *counter, int numsp)
1984
if(isTip(p->number, numsp))
1988
branches[*counter] = p->next;
1989
branches[*counter + 1] = p->next->next;
1991
*counter = *counter + 2;
1993
return ((2 + markBranches(branches, p->next->back, counter, numsp) +
1994
markBranches(branches, p->next->next->back, counter, numsp)));
2000
nodeptr findAnyTip(nodeptr p, int numsp)
2002
return isTip(p->number, numsp) ? p : findAnyTip(p->next->back, numsp);
2006
int randomInt(int n)
2011
void makePermutation(int *perm, int n, analdef *adef)
2017
for (i = 1; i <= n; i++)
2020
for (i = 1; i <= n; i++)
2022
k = (int)((double)(n + 1 - i) * randum(&adef->parsimonySeed));
2027
perm[i] = perm[i + k];
2039
boolean tipHomogeneityChecker(tree *tr, nodeptr p, int grouping)
2041
if(isTip(p->number, tr->mxtips))
2043
if(tr->constraintVector[p->number] != grouping)
2050
return (tipHomogeneityChecker(tr, p->next->back, grouping) && tipHomogeneityChecker(tr, p->next->next->back,grouping));
2060
void makeRandomTree(tree *tr, analdef *adef)
2062
nodeptr p, f, randomBranch;
2064
int *perm, branchCounter;
2067
branches = (nodeptr *)rax_malloc(sizeof(nodeptr) * (2 * tr->mxtips));
2068
perm = (int *)rax_malloc((tr->mxtips + 1) * sizeof(int));
2070
makePermutation(perm, tr->mxtips, adef);
2073
tr->nextnode = tr->mxtips + 1;
2075
buildSimpleTreeRandom(tr, perm[1], perm[2], perm[3]);
2077
while (tr->ntips < tr->mxtips)
2079
tr->bestParsimony = INT_MAX;
2080
nextsp = ++(tr->ntips);
2081
p = tr->nodep[perm[nextsp]];
2083
/*printf("ADDING SPECIES %d\n", nextsp);*/
2087
f = findAnyTip(tr->start, tr->mxtips);
2092
markBranches(branches, f, &branchCounter, tr->mxtips);
2094
assert(branchCounter == ((2 * (tr->ntips - 1)) - 3));
2096
randomBranch = branches[randomInt(branchCounter)];
2098
insertRandom(p->back, randomBranch, tr->numBranches);
2107
void nodeRectifier(tree *tr)
2109
nodeptr *np = (nodeptr *)rax_malloc(2 * tr->mxtips * sizeof(nodeptr));
2113
tr->start = tr->nodep[1];
2116
/* TODO why is tr->rooted set to FALSE here ?*/
2118
for(i = tr->mxtips + 1; i <= (tr->mxtips + tr->mxtips - 1); i++)
2119
np[i] = tr->nodep[i];
2121
reorderNodes(tr, np, tr->start->back, &count);
2128
void makeParsimonyTree(tree *tr, analdef *adef)
2130
makeParsimonyTreeFast(tr, adef, TRUE);
2133
void makeParsimonyTreeIncomplete(tree *tr, analdef *adef)
2135
makeParsimonyTreeFast(tr, adef, FALSE);
2139
static void setupBranchMetaInfo(tree *tr, nodeptr p, int nTips, branchInfo *bInf)
2142
countBranches = tr->branchCounter;
2144
if(isTip(p->number, tr->mxtips))
2146
p->bInf = &bInf[countBranches];
2147
p->back->bInf = &bInf[countBranches];
2149
bInf[countBranches].oP = p;
2150
bInf[countBranches].oQ = p->back;
2152
bInf[countBranches].epa->leftNodeNumber = p->number;
2153
bInf[countBranches].epa->rightNodeNumber = p->back->number;
2155
bInf[countBranches].epa->branchNumber = countBranches;
2156
bInf[countBranches].epa->originalBranchLength = p->z[0];
2158
tr->branchCounter = tr->branchCounter + 1;
2164
assert(p == p->next->next->next);
2166
p->bInf = &bInf[countBranches];
2167
p->back->bInf = &bInf[countBranches];
2169
bInf[countBranches].oP = p;
2170
bInf[countBranches].oQ = p->back;
2172
bInf[countBranches].epa->leftNodeNumber = p->number;
2173
bInf[countBranches].epa->rightNodeNumber = p->back->number;
2176
bInf[countBranches].epa->branchNumber = countBranches;
2177
bInf[countBranches].epa->originalBranchLength = p->z[0];
2179
tr->branchCounter = tr->branchCounter + 1;
2185
setupBranchMetaInfo(tr, q->back, nTips, bInf);
2195
static void setupJointFormat(tree *tr, nodeptr p, int ntips, branchInfo *bInf, int *count)
2197
if(isTip(p->number, tr->mxtips))
2199
p->bInf->epa->jointLabel = *count;
2200
*count = *count + 1;
2206
setupJointFormat(tr, p->next->back, ntips, bInf, count);
2207
setupJointFormat(tr, p->next->next->back, ntips, bInf, count);
2209
p->bInf->epa->jointLabel = *count;
2210
*count = *count + 1;
2221
static void setupBranchInfo(tree *tr, nodeptr q)
2224
originalNode = tr->nodep[tr->mxtips + 1];
2229
tr->branchCounter = 0;
2231
setupBranchMetaInfo(tr, q, tr->ntips, tr->bInf);
2233
assert(tr->branchCounter == tr->numberOfBranches);
2237
assert(tr->leftRootNode->back == tr->rightRootNode);
2238
assert(tr->leftRootNode == tr->rightRootNode->back);
2240
if(!isTip(tr->leftRootNode->number, tr->mxtips))
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);
2246
tr->leftRootNode->bInf->epa->jointLabel = count;
2247
tr->rootLabel = count;
2250
if(!isTip(tr->rightRootNode->number, tr->mxtips))
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);
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);
2263
assert(count == tr->numberOfBranches);
2266
static void testInsertFast(tree *tr, nodeptr r, nodeptr q)
2276
*inserts = tr->inserts;
2278
assert(!tr->grouped);
2280
hookupDefault(r->next, q, tr->numBranches);
2281
hookupDefault(r->next->next, x, tr->numBranches);
2283
newviewParsimony(tr, r);
2285
for(i = 0; i < tr->numberOfTipsForInsertion; i++)
2287
hookupDefault(r, tr->nodep[inserts[i]], tr->numBranches);
2289
tr->bestParsimony = INT_MAX;
2291
result = evaluateParsimony(tr, r, FALSE);
2293
r->back = (nodeptr) NULL;
2294
tr->nodep[inserts[i]]->back = (nodeptr) NULL;
2296
tr->bInf[q->bInf->epa->branchNumber].epa->parsimonyScore[i] = result;
2299
hookupDefault(q, x, tr->numBranches);
2301
r->next->next->back = r->next->back = (nodeptr) NULL;
2306
static void traverseTree(tree *tr, nodeptr r, nodeptr q)
2308
testInsertFast(tr, r, q);
2310
if(!isTip(q->number, tr->rdta->numsp))
2317
traverseTree(tr, r, a->back);
2327
unsigned int parsimonyScore;
2333
static int infoCompare(const void *p1, const void *p2)
2335
infoMP *rc1 = (infoMP *)p1;
2336
infoMP *rc2 = (infoMP *)p2;
2338
unsigned int i = rc1->parsimonyScore;
2339
unsigned int j = rc2->parsimonyScore;
2348
void classifyMP(tree *tr, analdef *adef)
2351
*informative = (int *)rax_malloc(sizeof(int) * (size_t)tr->cdta->endsite),
2364
jointFormatTreeFileName[1024],
2365
originalLabelledTreeFileName[1024],
2366
labelledTreeFileName[1024],
2367
likelihoodWeightsFileName[1024];
2370
*likelihoodWeightsFile,
2376
assert(adef->restart);
2378
determineUninformativeSites(tr, informative);
2380
compressDNA(tr, informative, TRUE);
2382
rax_free(informative);
2384
tr->ti = (int*)rax_malloc(sizeof(int) * 4 * (size_t)tr->mxtips);
2386
tr->numberOfBranches = 2 * tr->ntips - 3;
2388
printBothOpen("\nRAxML Evolutionary Placement Algorithm using parsimony\n");
2390
tr->bestParsimony = INT_MAX;
2392
score = evaluateParsimony(tr, tr->start->back, TRUE);
2394
printBothOpen("\nparsimony score of reference tree: %u\n\n", score);
2396
perm = (int *)rax_calloc(((size_t)tr->mxtips) + 1, sizeof(int));
2397
tr->inserts = (int *)rax_calloc((size_t)tr->mxtips, sizeof(int));
2399
markTips(tr->start, perm, tr->mxtips);
2400
markTips(tr->start->back, perm ,tr->mxtips);
2402
tr->numberOfTipsForInsertion = 0;
2404
for(i = 1; i <= tr->mxtips; i++)
2408
tr->inserts[tr->numberOfTipsForInsertion] = i;
2409
tr->numberOfTipsForInsertion = tr->numberOfTipsForInsertion + 1;
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);
2417
assert(tr->numberOfTipsForInsertion == (tr->mxtips - tr->ntips));
2419
tr->bInf = (branchInfo*)rax_malloc(tr->numberOfBranches * sizeof(branchInfo));
2421
for(i = 0; i < tr->numberOfBranches; i++)
2423
tr->bInf[i].epa = (epaBranchData*)rax_malloc(sizeof(epaBranchData));
2425
tr->bInf[i].epa->parsimonyScore = (unsigned int*)rax_malloc(tr->numberOfTipsForInsertion * sizeof(unsigned int));
2427
for(j = 0; j < tr->numberOfTipsForInsertion; j++)
2428
tr->bInf[i].epa->parsimonyScore[j] = INT_MAX;
2430
tr->bInf[i].epa->branchNumber = i;
2431
tr->bInf[i].epa->countThem = (int*)rax_calloc(tr->numberOfTipsForInsertion, sizeof(int));
2433
sprintf(tr->bInf[i].epa->branchLabel, "I%d", i);
2436
r = tr->nodep[(tr->nextnode)++];
2439
q = findAnyTip(tr->start, tr->rdta->numsp);
2441
assert(isTip(q->number, tr->rdta->numsp));
2442
assert(!isTip(q->back->number, tr->rdta->numsp));
2446
setupBranchInfo(tr, q);
2448
traverseTree(tr, r, q);
2450
printBothOpen("Overall Classification time: %f\n\n", gettime() - masterTime);
2453
strcpy(jointFormatTreeFileName, workdir);
2454
strcpy(originalLabelledTreeFileName, workdir);
2455
strcpy(labelledTreeFileName, workdir);
2456
strcpy(likelihoodWeightsFileName, workdir);
2458
strcat(jointFormatTreeFileName, "RAxML_portableTree.");
2459
strcat(originalLabelledTreeFileName, "RAxML_originalLabelledTree.");
2460
strcat(labelledTreeFileName, "RAxML_labelledTree.");
2461
strcat(likelihoodWeightsFileName, "RAxML_equallyParsimoniousPlacements.");
2463
strcat(jointFormatTreeFileName, run_id);
2464
strcat(originalLabelledTreeFileName, run_id);
2465
strcat(labelledTreeFileName, run_id);
2466
strcat(likelihoodWeightsFileName, run_id);
2468
strcat(jointFormatTreeFileName, ".jplace");
2470
rax_free(tr->tree_string);
2472
tr->treeStringLength *= 16;
2474
tr->tree_string = (char*)rax_calloc(tr->treeStringLength, sizeof(char));
2478
treeFile = myfopen(originalLabelledTreeFileName, "wb");
2479
Tree2StringClassify(tr->tree_string, tr, tr->inserts, TRUE, FALSE, FALSE);
2480
fprintf(treeFile, "%s\n", tr->tree_string);
2483
treeFile = myfopen(jointFormatTreeFileName, "wb");
2484
Tree2StringClassify(tr->tree_string, tr, tr->inserts, TRUE, TRUE, FALSE);
2486
fprintf(treeFile, "{\n");
2487
fprintf(treeFile, "\t\"tree\": \"%s\", \n", tr->tree_string);
2488
fprintf(treeFile, "\t\"placements\": [\n");
2491
inf = (infoMP*)rax_malloc(sizeof(infoMP) * tr->numberOfBranches);
2497
for(i = 0; i < tr->numberOfTipsForInsertion; i++)
2503
validEntries = tr->numberOfBranches;
2505
for(j = 0; j < tr->numberOfBranches; j++)
2507
inf[j].parsimonyScore = tr->bInf[j].epa->parsimonyScore[i];
2508
inf[j].number = tr->bInf[j].epa->jointLabel;
2511
qsort(inf, tr->numberOfBranches, sizeof(infoMP), infoCompare);
2515
lmax = inf[0].parsimonyScore;
2517
fprintf(treeFile, "\t{\"p\":[");
2519
while(j < validEntries && inf[j].parsimonyScore == lmax)
2523
if(tr->wasRooted && inf[j].number == tr->rootLabel)
2526
fprintf(treeFile, ",[%d, %u]", inf[j].number, inf[j].parsimonyScore);
2530
if(tr->wasRooted && inf[j].number == tr->rootLabel)
2533
fprintf(treeFile, "[%d, %u]", inf[j].number, inf[j].parsimonyScore);
2539
if(i == tr->numberOfTipsForInsertion - 1)
2540
fprintf(treeFile, "], \"n\":[\"%s\"]}\n", tr->nameList[tr->inserts[i]]);
2542
fprintf(treeFile, "], \"n\":[\"%s\"]},\n", tr->nameList[tr->inserts[i]]);
2545
fprintf(treeFile, "\t ],\n");
2546
/* fprintf(treeFile, "\t\"metadata\": {\"invocation\": \"RAxML EPA parsimony\"},\n");*/
2547
fprintf(treeFile, "\t\"metadata\": {\"invocation\": ");
2549
fprintf(treeFile, "\"");
2554
for(i = 0; i < globalArgc; i++)
2555
fprintf(treeFile,"%s ", globalArgv[i]);
2557
fprintf(treeFile, "\", \"raxml_version\": \"%s\"", programVersion);
2558
fprintf(treeFile,"},\n");
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");
2568
/* JSON format end */
2570
likelihoodWeightsFile = myfopen(likelihoodWeightsFileName, "wb");
2573
for(i = 0; i < tr->numberOfTipsForInsertion; i++)
2579
validEntries = tr->numberOfBranches;
2581
for(j = 0; j < tr->numberOfBranches; j++)
2583
inf[j].parsimonyScore = tr->bInf[j].epa->parsimonyScore[i];
2587
qsort(inf, tr->numberOfBranches, sizeof(infoMP), infoCompare);
2591
lmax = inf[0].parsimonyScore;
2593
while(j < validEntries && inf[j].parsimonyScore == lmax)
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;
2603
fclose(likelihoodWeightsFile);
2606
Tree2StringClassify(tr->tree_string, tr, tr->inserts, FALSE, FALSE, FALSE);
2607
treeFile = fopen(labelledTreeFileName, "wb");
2608
fprintf(treeFile, "%s\n", tr->tree_string);
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);
2622
#ifdef _SSE3_WAS_DEFINED
2626
#undef _SSE3_WAS_DEFINED
added
removed
GDE/RAxML/fastDNAparsimony.c
