source: trunk/GDE/RAxML/axml.h

Last change on this file was 19480, checked in by westram, 17 months ago
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 45.0 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
28 *  with thousands of taxa and mixed models".
29 *  Bioinformatics 2006; doi: 10.1093/bioinformatics/btl446
30 */
31
32
33#include <assert.h>
34#include <stdint.h>
35#include <stdio.h>
36
37#ifdef __AVX
38#define BYTE_ALIGNMENT 32
39#else
40#define BYTE_ALIGNMENT 16
41#endif
42
43
44#ifdef _USE_PTHREADS
45
46#include <pthread.h>
47
48#endif
49
50
51#define MAX_TIP_EV     0.999999999 /* max tip vector value, sum of EVs needs to be smaller than 1.0, otherwise the numerics break down */
52#define smoothings     32          /* maximum smoothing passes through tree */
53#define iterations     10          /* maximum iterations of iterations per insert */
54#define newzpercycle   1           /* iterations of makenewz per tree traversal */
55#define nmlngth        256         /* number of characters in species name */
56#define deltaz         0.00001     /* test of net branch length change in update */
57#define defaultz       0.9         /* value of z assigned as starting point */
58#define unlikely       -1.0E300    /* low likelihood for initialization */
59
60
61#define SUMMARIZE_LENGTH -3
62#define SUMMARIZE_LH     -2
63#define NO_BRANCHES      -1
64
65#define MASK_LENGTH 32
66#define VECTOR_LENGTH (NUM_BRANCHES / MASK_LENGTH)
67#define GET_BITVECTOR_LENGTH(x) ((x % MASK_LENGTH) ? (x / MASK_LENGTH + 1) : (x / MASK_LENGTH))
68
69#define zmin       1.0E-15  /* max branch prop. to -log(zmin) (= 34) */
70#define zmax (1.0 - 1.0E-6) /* min branch prop. to 1.0-zmax (= 1.0E-6) */
71
72#define twotothe256  \
73  115792089237316195423570985008687907853269984665640564039457584007913129639936.0
74                                                     /*  2**256 (exactly)  */
75
76#define minlikelihood  (1.0/twotothe256)
77#define minusminlikelihood -minlikelihood
78
79#define badRear         -1
80
81#define NUM_BRANCHES   128
82
83#define TRUE             1
84#define FALSE            0
85
86
87
88#define LIKELIHOOD_EPSILON 0.0000001
89
90#define THREAD_TO_DEBUG 1
91
92#define AA_SCALE 10.0
93#define AA_SCALE_PLUS_EPSILON 10.001
94
95/* ALPHA_MIN is critical -> numerical instability, eg for 4 discrete rate cats                    */
96/* and alpha = 0.01 the lowest rate r_0 is                                                        */
97/* 0.00000000000000000000000000000000000000000000000000000000000034878079110511010487             */
98/* which leads to numerical problems Table for alpha settings below:                              */
99/*                                                                                                */
100/* 0.010000 0.00000000000000000000000000000000000000000000000000000000000034878079110511010487    */
101/* 0.010000 yielded nasty numerical bugs in at least one case !                                   */
102/* 0.020000 0.00000000000000000000000000000044136090435925743185910935350715027016962154188875    */
103/* 0.030000 0.00000000000000000000476844846859006690412039180149775802624789852441798419292220    */
104/* 0.040000 0.00000000000000049522423236954066431210260930029681736928018820007024736185030633    */
105/* 0.050000 0.00000000000050625351310359203371872643495343928538368616365517027588794007897377    */
106/* 0.060000 0.00000000005134625283884191118711474021861409372524676086868566926568746566772461    */
107/* 0.070000 0.00000000139080650074206434685544624965062437960128249869740102440118789672851562    */
108/* 0.080000 0.00000001650681201563587066858709818343436959153791576682124286890029907226562500    */
109/* 0.090000 0.00000011301977332931251259273962858978301859735893231118097901344299316406250000    */
110/* 0.100000 0.00000052651925834844387815526344648331402709118265192955732345581054687500000000    */
111
112
113#define ALPHA_MIN    0.02
114#define ALPHA_MAX    1000.0
115
116#define RATE_MIN     0.0000001
117#define RATE_MAX     1000000.0
118
119#define INVAR_MIN    0.0001
120#define INVAR_MAX    0.9999
121
122#define TT_MIN       0.0000001
123#define TT_MAX       1000000.0
124
125#define FREQ_MIN     0.001
126
127#define LG4X_RATE_MIN 0.0000001
128#define LG4X_RATE_MAX 1000.0
129
130/*
131   previous values between 0.001 and 0.000001
132
133   TO AVOID NUMERICAL PROBLEMS WHEN FREQ == 0 IN PARTITIONED MODELS, ESPECIALLY WITH AA
134   previous value of FREQ_MIN was: 0.000001, but this seemed to cause problems with some
135   of the 7-state secondary structure models with some rather exotic small toy test datasets,
136   on the other hand 0.001 caused problems with some of the 16-state secondary structure models
137
138   For some reason the frequency settings seem to be repeatedly causing numerical problems
139   
140*/
141
142#define ITMAX 100
143
144
145
146#define SHFT(a,b,c,d)                (a)=(b);(b)=(c);(c)=(d);
147#define SIGN(a,b)                    ((b) > 0.0 ? fabs(a) : -fabs(a))
148
149#define ABS(x)    (((x)<0)   ?  (-(x)) : (x))
150#define MIN(x,y)  (((x)<(y)) ?    (x)  : (y))
151#define MAX(x,y)  (((x)>(y)) ?    (x)  : (y))
152#define NINT(x)   ((int) ((x)>0 ? ((x)+0.5) : ((x)-0.5)))
153
154
155#define LOG(x)  log(x)
156#define EXP(x)  exp(x)
157
158
159
160
161
162
163#define PointGamma(prob,alpha,beta)  PointChi2(prob,2.0*(alpha))/(2.0*(beta))
164
165#define programName        "RAxML"
166#define programVersion     "7.7.2"
167#define programDate        "July 31 2013"
168
169
170#define  TREE_EVALUATION                 0
171#define  BIG_RAPID_MODE                  1
172#define  CALC_BIPARTITIONS               2
173#define  SPLIT_MULTI_GENE                3
174#define  CHECK_ALIGNMENT                 4
175#define  PER_SITE_LL                     5
176#define  PARSIMONY_ADDITION              6
177#define  CLASSIFY_ML                     7
178#define  DISTANCE_MODE                   8
179#define  GENERATE_BS                     9
180#define  COMPUTE_ELW                     10
181#define  BOOTSTOP_ONLY                   11
182#define  COMPUTE_LHS                     12
183#define  COMPUTE_BIPARTITION_CORRELATION 13
184#define  COMPUTE_RF_DISTANCE             14
185#define  MORPH_CALIBRATOR                15
186#define  CONSENSUS_ONLY                  16
187#define  FAST_SEARCH                     17       
188#define  EPA_SITE_SPECIFIC_BIAS          18
189#define  SH_LIKE_SUPPORTS                19
190#define  CLASSIFY_MP                     20
191#define  ANCESTRAL_STATES                21
192#define  QUARTET_CALCULATION             22
193#define  THOROUGH_OPTIMIZATION           23
194#define  OPTIMIZE_BR_LEN_SCALER          24
195#define  ANCESTRAL_SEQUENCE_TEST         25
196#define  PLAUSIBILITY_CHECKER            26
197#define  CALC_BIPARTITIONS_IC            27
198#define  ROOT_TREE                       28
199
200#define M_GTRCAT         1
201#define M_GTRGAMMA       2
202#define M_BINCAT         3
203#define M_BINGAMMA       4
204#define M_PROTCAT        5
205#define M_PROTGAMMA      6
206#define M_32CAT          7
207#define M_32GAMMA        8
208#define M_64CAT          9
209#define M_64GAMMA        10
210
211#define DAYHOFF      0
212#define DCMUT        1
213#define JTT          2
214#define MTREV        3
215#define WAG          4
216#define RTREV        5
217#define CPREV        6
218#define VT           7
219#define BLOSUM62     8
220#define MTMAM        9
221#define LG           10
222#define MTART        11
223#define MTZOA        12
224#define PMB          13
225#define HIVB         14
226#define HIVW         15
227#define JTTDCMUT     16
228#define FLU          17
229#define DUMMY        18
230#define DUMMY2       19
231#define AUTO         20
232#define LG4          21
233#define LG4X         22
234#define PROT_FILE    23
235#define GTR_UNLINKED 24
236#define GTR          25  /* GTR always needs to be the last one */
237
238#define NUM_PROT_MODELS 26
239
240
241
242/* bipartition stuff */
243
244#define BIPARTITIONS_ALL       0
245#define GET_BIPARTITIONS_BEST  1
246#define DRAW_BIPARTITIONS_BEST 2
247#define BIPARTITIONS_BOOTSTOP  3
248#define BIPARTITIONS_RF  4
249#define GATHER_BIPARTITIONS_IC 5
250#define FIND_BIPARTITIONS_IC 6
251
252
253
254/* bootstopping stuff */
255
256#define BOOTSTOP_PERMUTATIONS 100
257#define START_BSTOP_TEST      10
258
259#define FC_THRESHOLD          99
260#define FC_SPACING            50
261#define FC_LOWER              0.99
262#define FC_INIT               20
263
264#define FREQUENCY_STOP 0
265#define MR_STOP        1
266#define MRE_STOP       2
267#define MRE_IGN_STOP   3
268
269#define MR_CONSENSUS 0
270#define MRE_CONSENSUS 1
271#define STRICT_CONSENSUS 2
272#define USER_DEFINED 3
273
274
275
276/* bootstopping stuff end */
277
278
279#define TIP_TIP     0
280#define TIP_INNER   1
281#define INNER_INNER 2
282
283#define MIN_MODEL        -1
284#define BINARY_DATA      0
285#define DNA_DATA         1
286#define AA_DATA          2
287#define SECONDARY_DATA   3
288#define SECONDARY_DATA_6 4
289#define SECONDARY_DATA_7 5
290#define GENERIC_32       6
291#define GENERIC_64       7
292#define MAX_MODEL        8
293
294#define SEC_6_A 0
295#define SEC_6_B 1
296#define SEC_6_C 2
297#define SEC_6_D 3
298#define SEC_6_E 4
299
300#define SEC_7_A 5
301#define SEC_7_B 6
302#define SEC_7_C 7
303#define SEC_7_D 8
304#define SEC_7_E 9
305#define SEC_7_F 10
306
307#define SEC_16   11
308#define SEC_16_A 12
309#define SEC_16_B 13
310#define SEC_16_C 14
311#define SEC_16_D 15
312#define SEC_16_E 16
313#define SEC_16_F 17
314#define SEC_16_I 18
315#define SEC_16_J 19
316#define SEC_16_K 20
317
318#define ORDERED_MULTI_STATE 0
319#define MK_MULTI_STATE      1
320#define GTR_MULTI_STATE     2
321
322
323
324
325
326#define CAT         0
327#define GAMMA       1
328#define GAMMA_I     2
329
330
331
332
333typedef  int boolean;
334
335
336typedef struct {
337  double lh;
338  int tree;
339  double weight;
340} elw;
341
342struct ent
343{
344  unsigned int *bitVector;
345  unsigned int *treeVector;
346  unsigned int amountTips;
347  int *supportVector;
348  unsigned int bipNumber;
349  unsigned int bipNumber2;
350  unsigned int supportFromTreeset[2]; 
351  struct ent *next;
352};
353
354typedef struct ent entry;
355
356typedef unsigned int hashNumberType;
357
358typedef unsigned int parsimonyNumber;
359
360/*typedef uint_fast32_t parsimonyNumber;*/
361
362#define PCF 32
363
364/*
365  typedef uint64_t parsimonyNumber;
366
367  #define PCF 16
368
369
370typedef unsigned char parsimonyNumber;
371
372#define PCF 2
373*/
374
375typedef struct
376{
377  hashNumberType tableSize;
378  entry **table;
379  hashNumberType entryCount;
380}
381  hashtable;
382
383
384struct stringEnt
385{
386  int nodeNumber;
387  char *word;
388  struct stringEnt *next;
389};
390
391typedef struct stringEnt stringEntry;
392 
393typedef struct
394{
395  hashNumberType tableSize;
396  stringEntry **table;
397}
398  stringHashtable;
399
400
401
402
403typedef struct ratec
404{
405  double accumulatedSiteLikelihood;
406  double rate;
407}
408  rateCategorize;
409
410
411typedef struct
412{
413  int tipCase;
414  int pNumber;
415  int qNumber;
416  int rNumber;
417  double qz[NUM_BRANCHES];
418  double rz[NUM_BRANCHES];
419} traversalInfo;
420
421typedef struct
422{
423  traversalInfo *ti;
424  int count;
425} traversalData;
426
427
428struct noderec;
429
430typedef struct epBrData
431{
432  int    *countThem;
433  int    *executeThem;
434  unsigned int *parsimonyScore;
435  double *branches;
436  double *distalBranches; 
437  double *likelihoods;
438  double originalBranchLength;
439  char branchLabel[64];
440  int leftNodeNumber;
441  int rightNodeNumber;
442  int *leftScaling;
443  int *rightScaling;
444  double branchLengths[NUM_BRANCHES];
445  double *left;
446  double *right;
447  int branchNumber;
448  int jointLabel;
449} epaBranchData;
450
451typedef struct
452{
453  epaBranchData *epa;
454
455  unsigned int *vector; 
456  int support;
457  int *supports;
458  double ic;
459  double icAll;
460  struct noderec *oP;
461  struct noderec *oQ;
462} branchInfo;
463
464
465
466
467
468
469
470
471typedef struct
472{
473  boolean valid;
474  int partitions;
475  int *partitionList;
476}
477  linkageData;
478
479typedef struct
480{
481  int entries;
482  linkageData* ld;
483}
484  linkageList;
485
486
487typedef  struct noderec
488{ 
489  branchInfo      *bInf;
490  double           z[NUM_BRANCHES];
491  struct noderec  *next;
492  struct noderec  *back;
493  hashNumberType   hash;
494  int              support;
495  int              number;
496  char             x;
497}
498  node, *nodeptr;
499
500typedef struct
501  {
502    double lh;
503    double pendantBranch;
504    double distalBranch;   
505    int number;
506  }
507  info;
508
509typedef struct bInf {
510  double likelihood;
511  nodeptr node;
512} bestInfo;
513
514typedef struct iL {
515  bestInfo *list;
516  int n;
517  int valid;
518} infoList;
519
520
521
522
523typedef  struct
524{
525  int              numsp;
526  int              sites;
527  unsigned char             **y;
528  unsigned char             *y0;
529  unsigned char             *yBUF;
530  int              *wgt;
531} rawdata;
532
533typedef  struct {
534  int             *alias;       /* site representing a pattern */
535  int             *aliaswgt;    /* weight by pattern */
536  int             *rateCategory;
537  int              endsite;     /* # of sequence patterns */
538  double          *patrat;      /* rates per pattern */
539  double          *patratStored; 
540} cruncheddata;
541
542
543
544
545typedef struct {
546  int     states;
547  int     maxTipStates;
548  size_t    lower;
549  size_t     upper;
550  size_t     width;
551  int     dataType;
552  int     protModels;
553  int     autoProtModels;
554  boolean usePredefinedProtFreqs;
555  int     mxtips;
556  int     numberOfCategories;
557  int             **expVector;
558  double          **xVector;
559  size_t             *xSpaceVector;
560  size_t             *expSpaceVector;
561 
562  unsigned char            **yVector;
563 
564  char   *partitionName;
565  char   proteinSubstitutionFileName[2048];
566  double externalAAMatrix[420];
567
568  double *sumBuffer;
569   double *gammaRates;
570
571  double *EIGN;
572  double *EV;
573  double *EI; 
574
575  double *left;
576  double *right;
577
578  /* LG4 */
579
580  double *EIGN_LG4[4];
581  double *EV_LG4[4];
582  double *EI_LG4[4]; 
583
584 
585
586  double *frequencies_LG4[4];
587  double *tipVector_LG4[4];
588  double *substRates_LG4[4];
589 
590  /* LG4X */
591
592  double weights[4];
593  double weightExponents[4];
594
595  /* LG4 */
596
597  double *frequencies;
598  double *tipVector;
599 
600  double *substRates;
601  double *perSiteLL;
602 
603  double *perSiteRates;
604  double *unscaled_perSiteRates;
605
606  unsigned int    *globalScaler;
607 
608  int    *wgt;
609  int    *invariant;
610  int    *rateCategory;
611  int    *symmetryVector;
612  int    *frequencyGrouping;
613  boolean nonGTR;
614  double alpha;
615  double propInvariant;
616
617  int gapVectorLength;
618  unsigned int *gapVector;
619  double *gapColumn;
620
621  size_t initialGapVectorSize;
622
623  size_t parsimonyLength;
624  parsimonyNumber *parsVect; 
625
626  double brLenScaler;
627
628} pInfo;
629
630
631
632typedef struct 
633{
634  int left;
635  int right;
636  double likelihood;
637} lhEntry;
638
639
640typedef struct 
641{
642  int count;
643  int size;
644  lhEntry *entries;
645} lhList;
646
647
648typedef struct List_{
649  void *value;                 
650  struct List_ *next; 
651} List;
652
653
654
655
656typedef  struct  {
657  boolean optimizeAllTrees;
658
659 
660  boolean saveMemory;
661 
662  int    *resample;
663
664  int numberOfBranches;
665  int    numberOfTipsForInsertion;
666  int    *readPartition;
667  boolean perPartitionEPA;
668  int    *inserts;
669  int    branchCounter;
670 
671  int *ti; 
672 
673  int numberOfTrees; 
674
675  stringHashtable  *nameHash;
676
677  pInfo            *partitionData;
678  pInfo            *initialPartitionData;
679  pInfo            *extendedPartitionData;
680
681  int              *dataVector;
682  int              *initialDataVector;
683  int              *extendedDataVector;
684
685  int              *patternPosition;
686  int              *columnPosition;
687
688  char             *secondaryStructureInput;
689
690  boolean          *executeModel;
691
692  double           *perPartitionLH;
693  double           *storedPerPartitionLH;
694
695  traversalData td[1];
696
697  unsigned int *parsimonyScore;
698
699  int              maxCategories;
700
701  double           *sumBuffer;
702  double           *perSiteLL; 
703  double           coreLZ[NUM_BRANCHES];
704  int              modelNumber;
705  int              multiBranch;
706  int              numBranches;
707  int              bootStopCriterion;
708  int              consensusType;
709  int              consensusUserThreshold;
710  double           wcThreshold;
711
712 
713  double          *storedBrLens;
714  boolean         useBrLenScaler;
715 
716
717 
718 
719 
720
721  boolean          useFastScaling;
722 
723  branchInfo       *bInf;
724
725  int              multiStateModel;
726
727
728  size_t innerNodes;
729
730  boolean curvatOK[NUM_BRANCHES];
731  /* the stuff below is shared among DNA and AA, span does
732     not change depending on datatype */
733
734  double           *invariants;
735  double           *fracchanges;
736 
737  double           *rawFracchanges;
738
739  /* model stuff end */
740
741  unsigned char             **yVector;
742  int              secondaryStructureModel;
743  int              discreteRateCategories;
744  int              originalCrunchedLength;
745  int              fullSites;
746  int              *originalModel;
747  int              *originalDataVector;
748  int              *originalWeights;
749  int              *secondaryStructurePairs;
750
751
752  double            *partitionContributions;
753
754  double            fracchange;
755  double            rawFracchange;
756
757  double            lhCutoff;
758  double            lhAVG;
759  unsigned long     lhDEC;
760  unsigned long     itCount;
761  int               numberOfInvariableColumns;
762  int               weightOfInvariableColumns;
763  int               rateHetModel;
764
765  double           startLH;
766  double           endLH;
767  double           likelihood;
768  double          *likelihoods;
769  int             *invariant;
770  node           **nodep;
771  node            *start;
772  int              mxtips;
773  int              mxtipsVector[NUM_BRANCHES];
774  int              *model;
775
776  int              *constraintVector;
777  int              numberOfSecondaryColumns;
778  boolean          searchConvergenceCriterion;
779  int              branchLabelCounter;
780  int              ntips;
781  int              nextnode;
782  int              NumberOfModels;
783  int              parsimonyLength;
784 
785  int              checkPointCounter;
786  int              treeID;
787  int              numberOfOutgroups;
788  int             *outgroupNums;
789  char           **outgroups;
790  boolean          useEpaHeuristics;
791  double           fastEPAthreshold;
792  boolean          bigCutoff;
793  boolean          partitionSmoothed[NUM_BRANCHES];
794  boolean          partitionConverged[NUM_BRANCHES];
795  boolean          rooted;
796  boolean          grouped;
797  boolean          constrained;
798  boolean          doCutoff;
799  boolean          catOnly;
800  rawdata         *rdta;
801  cruncheddata    *cdta;
802
803  char **nameList;
804  char *tree_string;
805  int treeStringLength;
806  unsigned int bestParsimony;
807  double bestOfNode;
808  nodeptr removeNode;
809  nodeptr insertNode;
810
811  double zqr[NUM_BRANCHES];
812  double currentZQR[NUM_BRANCHES];
813
814  double currentLZR[NUM_BRANCHES];
815  double currentLZQ[NUM_BRANCHES];
816  double currentLZS[NUM_BRANCHES];
817  double currentLZI[NUM_BRANCHES];
818  double lzs[NUM_BRANCHES];
819  double lzq[NUM_BRANCHES];
820  double lzr[NUM_BRANCHES];
821  double lzi[NUM_BRANCHES];
822
823 
824  int mr_thresh;
825
826  boolean wasRooted;
827  nodeptr leftRootNode;
828  nodeptr rightRootNode;
829  int rootLabel;
830 
831  boolean useGammaMedian;
832  boolean noRateHet;
833
834#ifdef _USE_PTHREADS
835
836  double *ancestralStates;
837
838  hashtable *h;
839 
840   
841   
842
843  double *temporaryVector; 
844  int    *temporaryScaling;
845  double *temporarySumBuffer;
846  size_t contiguousVectorLength;
847  size_t contiguousScalingLength; 
848 
849 
850
851  int *contiguousRateCategory;
852  int *contiguousWgt;
853  int *contiguousInvariant; 
854
855  unsigned char **contiguousTips;
856 
857 
858  unsigned char *y_ptr;
859
860 
861
862  double *perSiteLLPtr;
863  int    *wgtPtr;
864  int    *invariantPtr;
865  int    *rateCategoryPtr;
866
867  int threadID;
868  double lower_spacing;
869  double upper_spacing;
870  double *lhs;
871
872  /* stuff for parallel MRE */
873
874  /* added by aberer */ 
875  entry **sbi;
876  entry **sbw;
877  int *len;   
878
879  /* mre */
880  int sectionEnd;
881  int bipStatusLen;
882  entry **entriesOfSection;
883  int recommendedAmountJobs;
884 
885
886  /* used for printBip */
887  boolean *hasAncestor; 
888  List **listOfDirectChildren;
889  unsigned int bitVectorLength;                 /* we now need this also in sequential mode */
890  entry **consensusBips;
891  int consensusBipLen;          /* also used in mre */ 
892  int maxBips;
893  int *bipStatus;
894
895
896  /* for parallel hash insert */
897  pthread_mutex_t** mutexesForHashing; 
898
899
900 
901
902#endif
903
904
905  int *origNumSitePerModel;
906
907} tree;
908
909
910/***************************************************************/
911
912typedef struct
913{
914  double z[NUM_BRANCHES];
915  nodeptr p, q;
916  int cp, cq;
917}
918  connectRELL, *connptrRELL;
919
920typedef  struct
921{
922  connectRELL     *connect; 
923  int             start;
924  double          likelihood;
925}
926  topolRELL;
927
928
929typedef  struct
930{
931  int max;
932  topolRELL **t;
933}
934  topolRELL_LIST;
935
936
937/**************************************************************/
938
939
940
941typedef struct conntyp {
942    double           z[NUM_BRANCHES];           /* branch length */
943    node            *p, *q;       /* parent and child sectors */
944    void            *valptr;      /* pointer to value of subtree */
945    int              descend;     /* pointer to first connect of child */
946    int              sibling;     /* next connect from same parent */
947    } connect, *connptr;
948
949typedef  struct {
950    double           likelihood;
951  int              initialTreeNumber;
952    connect         *links;       /* pointer to first connect (start) */
953    node            *start;
954    int              nextlink;    /* index of next available connect */
955                                  /* tr->start = tpl->links->p */
956    int              ntips;
957    int              nextnode;
958    int              scrNum;      /* position in sorted list of scores */
959    int              tplNum;      /* position in sorted list of trees */
960
961    } topol;
962
963typedef struct {
964    double           best;        /* highest score saved */
965    double           worst;       /* lowest score saved */
966    topol           *start;       /* starting tree for optimization */
967    topol          **byScore;
968    topol          **byTopol;
969    int              nkeep;       /* maximum topologies to save */
970    int              nvalid;      /* number of topologies saved */
971    int              ninit;       /* number of topologies initialized */
972    int              numtrees;    /* number of alternatives tested */
973    boolean          improved;
974    } bestlist;
975
976#define PHYLIP 0
977#define FASTA  1
978
979typedef  struct {
980  int              categories;
981  int              model;
982  int              bestTrav;
983  int              max_rearrange;
984  int              stepwidth;
985  int              initial;
986  boolean          initialSet;
987  int              mode;
988  long             boot;
989  long             rapidBoot;
990  boolean          bootstrapBranchLengths;
991  boolean          restart;
992  boolean          useWeightFile;
993  boolean          useMultipleModel;
994  boolean          constraint;
995  boolean          grouping;
996  boolean          randomStartingTree;
997  boolean          useInvariant;
998  int            protEmpiricalFreqs;
999  int            proteinMatrix;
1000  int            checkpoints;
1001  int            startingTreeOnly;
1002  int            multipleRuns;
1003  long           parsimonySeed;
1004  boolean        perGeneBranchLengths;
1005  boolean        likelihoodTest;
1006  boolean        outgroup;
1007  boolean        permuteTreeoptimize;
1008  boolean        allInOne;
1009  boolean        generateBS;
1010  boolean        bootStopping;
1011  boolean        useExcludeFile;
1012  boolean        userProteinModel;
1013  boolean        computeELW;
1014  boolean        computeDistance;
1015  boolean        compressPatterns;
1016  boolean        useSecondaryStructure; 
1017  double         likelihoodEpsilon;
1018  double         gapyness;
1019  int            similarityFilterMode;
1020  double        externalAAMatrix[420];
1021  boolean       useFloat;
1022  boolean       readTaxaOnly;
1023  boolean       veryFast;
1024  boolean       useBinaryModelFile;
1025  boolean       leaveDropMode;
1026  int           slidingWindowSize;
1027  boolean       checkForUndeterminedSequences;
1028  boolean       useQuartetGrouping;
1029  int           alignmentFileType;
1030  boolean       calculateIC;
1031  boolean       verboseIC;
1032  boolean       stepwiseAdditionOnly;
1033} analdef;
1034
1035
1036typedef struct 
1037{
1038  int leftLength;
1039  int rightLength;
1040  int eignLength;
1041  int evLength;
1042  int eiLength;
1043  int substRatesLength;
1044  int frequenciesLength;
1045  int tipVectorLength;
1046  int symmetryVectorLength;
1047  int frequencyGroupingLength;
1048
1049  boolean nonGTR;
1050
1051  int undetermined;
1052
1053  const char *inverseMeaning;
1054
1055  int states;
1056
1057  boolean smoothFrequencies;
1058
1059  const unsigned  int *bitVector;
1060
1061} partitionLengths;
1062
1063/****************************** FUNCTIONS ****************************************************/
1064
1065extern void computePlacementBias(tree *tr, analdef *adef);
1066
1067extern int lookupWord(char *s, stringHashtable *h);
1068
1069extern void getDataTypeString(tree *tr, int model, char typeOfData[1024]);
1070
1071extern unsigned int genericBitCount(unsigned int* bitVector, unsigned int bitVectorLength);
1072extern int countTips(nodeptr p, int numsp);
1073extern entry *initEntry(void);
1074extern void computeRogueTaxa(tree *tr, char* treeSetFileName, analdef *adef);
1075extern unsigned int precomputed16_bitcount(unsigned int n);
1076
1077#define BIT_COUNT(x)  precomputed16_bitcount(x)
1078
1079
1080
1081
1082
1083
1084extern partitionLengths * getPartitionLengths(pInfo *p);
1085extern boolean getSmoothFreqs(int dataType);
1086extern const unsigned int *getBitVector(int dataType);
1087extern unsigned char getUndetermined(int dataType);
1088extern int getStates(int dataType);
1089extern char getInverseMeaning(int dataType, unsigned char state);
1090extern void printModelParams(tree *tr, analdef *adef);
1091extern double gettime ( void );
1092extern int gettimeSrand ( void );
1093extern double randum ( long *seed );
1094
1095extern void getxnode ( nodeptr p );
1096extern void hookup ( nodeptr p, nodeptr q, double *z, int numBranches);
1097extern void hookupDefault ( nodeptr p, nodeptr q, int numBranches);
1098extern boolean whitechar ( int ch );
1099extern void errorExit ( int e ) __attribute__((noreturn));
1100extern void printResult ( tree *tr, analdef *adef, boolean finalPrint );
1101extern void printBootstrapResult ( tree *tr, analdef *adef, boolean finalPrint );
1102extern void printBipartitionResult ( tree *tr, analdef *adef, boolean finalPrint, boolean printIC);
1103extern void printLog ( tree *tr, analdef *adef, boolean finalPrint );
1104extern void printStartingTree ( tree *tr, analdef *adef, boolean finalPrint );
1105extern void writeInfoFile ( analdef *adef, tree *tr, double t );
1106extern int main ( int argc, char *argv[] );
1107extern void calcBipartitions ( tree *tr, analdef *adef, char *bestTreeFileName, char *bootStrapFileName );
1108extern void calcBipartitions_IC ( tree *tr, analdef *adef, char *bestTreeFileName, char *bootStrapFileName );
1109
1110extern void initReversibleGTR (tree *tr, int model);
1111extern double LnGamma ( double alpha );
1112extern double IncompleteGamma ( double x, double alpha, double ln_gamma_alpha );
1113extern double PointNormal ( double prob );
1114extern double PointChi2 ( double prob, double v );
1115extern void makeGammaCats (double alpha, double *gammaRates, int K,  boolean useMedian);
1116extern void initModel ( tree *tr, rawdata *rdta, cruncheddata *cdta, analdef *adef );
1117extern void doAllInOne ( tree *tr, analdef *adef );
1118
1119extern void classifyML(tree *tr, analdef *adef);
1120extern void classifyMP(tree *tr, analdef *adef);
1121extern void markTips(nodeptr p, int *perm, int maxTips);
1122extern char *Tree2StringClassify(char *treestr, tree *tr, int *inserts, 
1123                                 boolean  originalTree, boolean jointLabels, boolean likelihood);
1124
1125
1126extern void doBootstrap ( tree *tr, analdef *adef, rawdata *rdta, cruncheddata *cdta );
1127extern void doInference ( tree *tr, analdef *adef, rawdata *rdta, cruncheddata *cdta );
1128extern void resetBranches ( tree *tr );
1129extern void scaleBranches(tree *tr, boolean fromFile);
1130extern void modOpt ( tree *tr, analdef *adef , boolean resetModel, double likelihoodEpsilon);
1131
1132extern void plausibilityChecker(tree *tr, analdef *adef);
1133
1134extern void parsePartitions ( analdef *adef, rawdata *rdta, tree *tr);
1135extern void computeBOOTRAPID (tree *tr, analdef *adef, long *radiusSeed);
1136extern void optimizeRAPID ( tree *tr, analdef *adef );
1137extern void thoroughOptimization ( tree *tr, analdef *adef, topolRELL_LIST *rl, int index );
1138extern int treeOptimizeThorough ( tree *tr, int mintrav, int maxtrav);
1139
1140extern int checker ( tree *tr, nodeptr p );
1141extern int randomInt ( int n );
1142extern void makePermutation ( int *perm, int n, analdef *adef );
1143extern boolean tipHomogeneityChecker ( tree *tr, nodeptr p, int grouping );
1144extern void makeRandomTree ( tree *tr, analdef *adef );
1145extern void nodeRectifier ( tree *tr );
1146extern void makeParsimonyTreeThorough(tree *tr, analdef *adef);
1147extern void makeParsimonyTree ( tree *tr, analdef *adef );
1148extern void makeParsimonyTreeFast(tree *tr, analdef *adef, boolean full);
1149extern void makeParsimonyTreeIncomplete ( tree *tr, analdef *adef );
1150extern void makeParsimonyInsertions(tree *tr, nodeptr startNodeQ, nodeptr startNodeR);
1151
1152
1153
1154extern FILE *myfopen(const char *path, const char *mode);
1155
1156
1157extern boolean initrav ( tree *tr, nodeptr p );
1158extern void initravPartition ( tree *tr, nodeptr p, int model );
1159extern boolean update ( tree *tr, nodeptr p );
1160extern boolean smooth ( tree *tr, nodeptr p );
1161extern boolean smoothTree ( tree *tr, int maxtimes );
1162extern boolean localSmooth ( tree *tr, nodeptr p, int maxtimes );
1163extern boolean localSmoothMulti(tree *tr, nodeptr p, int maxtimes, int model);
1164extern void initInfoList ( int n );
1165extern void freeInfoList ( void );
1166extern void insertInfoList ( nodeptr node, double likelihood );
1167extern boolean smoothRegion ( tree *tr, nodeptr p, int region );
1168extern boolean regionalSmooth ( tree *tr, nodeptr p, int maxtimes, int region );
1169extern nodeptr removeNodeBIG ( tree *tr, nodeptr p, int numBranches);
1170extern nodeptr removeNodeRestoreBIG ( tree *tr, nodeptr p );
1171extern boolean insertBIG ( tree *tr, nodeptr p, nodeptr q, int numBranches);
1172extern boolean insertRestoreBIG ( tree *tr, nodeptr p, nodeptr q );
1173extern boolean testInsertBIG ( tree *tr, nodeptr p, nodeptr q );
1174extern void addTraverseBIG ( tree *tr, nodeptr p, nodeptr q, int mintrav, int maxtrav );
1175extern int rearrangeBIG ( tree *tr, nodeptr p, int mintrav, int maxtrav );
1176extern void traversalOrder ( nodeptr p, int *count, nodeptr *nodeArray );
1177extern double treeOptimizeRapid ( tree *tr, int mintrav, int maxtrav, analdef *adef, bestlist *bt);
1178extern boolean testInsertRestoreBIG ( tree *tr, nodeptr p, nodeptr q );
1179extern void restoreTreeFast ( tree *tr );
1180extern int determineRearrangementSetting ( tree *tr, analdef *adef, bestlist *bestT, bestlist *bt );
1181extern void computeBIGRAPID ( tree *tr, analdef *adef, boolean estimateModel);
1182extern boolean treeEvaluate ( tree *tr, double smoothFactor );
1183extern boolean treeEvaluatePartition ( tree *tr, double smoothFactor, int model );
1184
1185extern void initTL ( topolRELL_LIST *rl, tree *tr, int n );
1186extern void freeTL ( topolRELL_LIST *rl);
1187extern void restoreTL ( topolRELL_LIST *rl, tree *tr, int n );
1188extern void resetTL ( topolRELL_LIST *rl );
1189extern void saveTL ( topolRELL_LIST *rl, tree *tr, int index );
1190
1191extern int  saveBestTree (bestlist *bt, tree *tr);
1192extern int  recallBestTree (bestlist *bt, int rank, tree *tr);
1193extern int initBestTree ( bestlist *bt, int newkeep, int numsp );
1194extern void resetBestTree ( bestlist *bt );
1195extern boolean freeBestTree ( bestlist *bt );
1196
1197
1198extern char *Tree2String ( char *treestr, tree *tr, nodeptr p, boolean printBranchLengths, boolean printNames, boolean printLikelihood, 
1199                           boolean rellTree, boolean finalPrint, analdef *adef, int perGene, boolean branchLabelSupport, boolean printSHSupport, boolean printIC, boolean printSHSupports);
1200extern void printTreePerGene(tree *tr, analdef *adef, const char *fileName, const char *permission);
1201
1202
1203extern int treeFindTipName(FILE *fp, tree *tr, boolean check);
1204extern int treeReadLen (FILE *fp, tree *tr, boolean readBranches, boolean readNodeLabels, boolean topologyOnly, analdef *adef, boolean completeTree, boolean storeBranchLabels);
1205extern boolean treeReadLenMULT ( FILE *fp, tree *tr, analdef *adef );
1206
1207extern int readMultifurcatingTree(FILE *fp, tree *tr, analdef *adef);
1208extern void freeMultifurcations(tree *tr);
1209extern void allocateMultifurcations(tree *tr, tree *smallTree);
1210
1211extern void getStartingTree ( tree *tr, analdef *adef );
1212extern double treeLength(tree *tr, int model);
1213
1214extern void computeBootStopOnly(tree *tr, char *bootStrapFileName, analdef *adef);
1215extern boolean bootStop(tree *tr, hashtable *h, int numberOfTrees, double *pearsonAverage, unsigned int **bitVectors, int treeVectorLength, unsigned int vectorLength);
1216extern void computeConsensusOnly(tree *tr, char* treeSetFileName, analdef *adef, boolean computeIC);
1217extern double evaluatePartialGeneric (tree *, int i, double ki, int _model);
1218extern double evaluateGeneric (tree *tr, nodeptr p);
1219extern void newviewGeneric (tree *tr, nodeptr p);
1220extern void newviewGenericMulti (tree *tr, nodeptr p, int model);
1221extern void newviewGenericMasked (tree *tr, nodeptr p);
1222extern void makenewzGeneric(tree *tr, nodeptr p, nodeptr q, double *z0, int maxiter, double *result, boolean mask);
1223extern void makenewzGenericDistance(tree *tr, int maxiter, double *z0, double *result, int taxon1, int taxon2);
1224extern double evaluatePartitionGeneric (tree *tr, nodeptr p, int model);
1225extern void newviewPartitionGeneric (tree *tr, nodeptr p, int model);
1226extern double evaluateGenericVector (tree *tr, nodeptr p);
1227extern void categorizeGeneric (tree *tr, nodeptr p);
1228extern double makenewzPartitionGeneric(tree *tr, nodeptr p, nodeptr q, double z0, int maxiter, int model);
1229extern boolean isTip(int number, int maxTips);
1230extern void computeTraversalInfo(nodeptr p, traversalInfo *ti, int *counter, int maxTips, int numBranches);
1231
1232
1233
1234extern void   newviewIterative(tree *);
1235
1236extern double evaluateIterative(tree *, boolean writeVector);
1237
1238extern double FABS(double x);
1239
1240
1241
1242
1243extern void makenewzIterative(tree *);
1244extern void execCore(tree *, volatile double *dlnLdlz, volatile double *d2lnLdlz2);
1245
1246
1247
1248extern void determineFullTraversal(nodeptr p, tree *tr);
1249/*extern void optRateCat(tree *, int i, double lower_spacing, double upper_spacing, double *lhs);*/
1250
1251extern unsigned int evaluateParsimonyIterative(tree *);
1252extern void newviewParsimonyIterative(tree *);
1253
1254extern unsigned int evaluateParsimonyIterativeFast(tree *);
1255extern void newviewParsimonyIterativeFast(tree *);
1256
1257
1258extern void initravParsimonyNormal(tree *tr, nodeptr p);
1259
1260extern double evaluateGenericInitrav (tree *tr, nodeptr p);
1261extern double evaluateGenericInitravPartition(tree *tr, nodeptr p, int model);
1262extern void evaluateGenericVectorIterative(tree *, int startIndex, int endIndex);
1263extern void categorizeIterative(tree *, int startIndex, int endIndex);
1264extern void onlyInitrav(tree *tr, nodeptr p);
1265extern void onlyInitravPartition(tree *tr, nodeptr p, int model);
1266
1267extern void fixModelIndices(tree *tr, int endsite, boolean fixRates);
1268extern void calculateModelOffsets(tree *tr);
1269extern void gammaToCat(tree *tr);
1270extern void catToGamma(tree *tr, analdef *adef);
1271extern void handleExcludeFile(tree *tr, analdef *adef, rawdata *rdta);
1272extern void printBaseFrequencies(tree *tr);
1273extern nodeptr findAnyTip(nodeptr p, int numsp);
1274
1275extern void parseProteinModel(double *externalAAMatrix, char *fileName);
1276extern int filexists(char *filename);
1277extern void computeFullTraversalInfo(nodeptr p, traversalInfo *ti, int *counter, int maxTips, int numBranches);
1278
1279extern void computeNextReplicate(tree *tr, long *seed, int *originalRateCategories, int *originalInvariant, boolean isRapid, boolean fixRates);
1280/*extern void computeNextReplicate(tree *tr, analdef *adef, int *originalRateCategories, int *originalInvariant);*/
1281
1282extern void reductionCleanup(tree *tr, int *originalRateCategories, int *originalInvariant);
1283extern void parseSecondaryStructure(tree *tr, analdef *adef, int sites);
1284extern void printPartitions(tree *tr);
1285extern void calcDiagptable(double z, int data, int numberOfCategories, double *rptr, double *EIGN, double *diagptable);
1286extern void compareBips(tree *tr, char *bootStrapFileName, analdef *adef);
1287extern void computeRF(tree *tr, char *bootStrapFileName, analdef *adef);
1288
1289
1290extern  unsigned int **initBitVector(tree *tr, unsigned int *vectorLength);
1291extern hashtable *copyHashTable(hashtable *src, unsigned int vectorLength);
1292extern hashtable *initHashTable(unsigned int n);
1293extern void cleanupHashTable(hashtable *h, int state);
1294extern double convergenceCriterion(hashtable *h, int mxtips);
1295extern void freeBitVectors(unsigned int **v, int n);
1296extern void freeHashTable(hashtable *h);
1297extern stringHashtable *initStringHashTable(hashNumberType n);
1298extern void addword(char *s, stringHashtable *h, int nodeNumber);
1299
1300
1301extern void printBothOpen(const char* format, ... );
1302extern void printBothOpenMPI(const char* format, ... );
1303extern void initRateMatrix(tree *tr);
1304
1305extern void bitVectorInitravSpecial(unsigned int **bitVectors, nodeptr p, int numsp, unsigned int vectorLength, hashtable *h, int treeNumber, int function, branchInfo *bInf,
1306                                    int *countBranches, int treeVectorLength, boolean traverseOnly, boolean computeWRF);
1307
1308extern int getIncrement(tree *tr, int model);
1309
1310extern void fastSearch(tree *tr, analdef *adef, rawdata *rdta, cruncheddata *cdta);
1311extern void shSupports(tree *tr, analdef *adef, rawdata *rdta, cruncheddata *cdta);
1312
1313extern FILE *getNumberOfTrees(tree *tr, char *fileName, analdef *adef);
1314
1315extern void writeBinaryModel(tree *tr);
1316extern void readBinaryModel(tree *tr);
1317extern void treeEvaluateRandom (tree *tr, double smoothFactor);
1318extern void treeEvaluateProgressive(tree *tr);
1319
1320
1321
1322extern boolean issubset(unsigned int* bipA, unsigned int* bipB, unsigned int vectorLen, unsigned int firstIndex);
1323extern boolean compatible(entry* e1, entry* e2, unsigned int bvlen);
1324
1325extern void nniSmooth(tree *tr, nodeptr p, int maxtimes);
1326
1327extern int *permutationSH(tree *tr, int nBootstrap, long _randomSeed);
1328
1329extern void updatePerSiteRates(tree *tr, boolean scaleRates);
1330
1331extern void newviewIterativeAncestral(tree *tr);
1332extern void newviewGenericAncestral(tree *tr, nodeptr p, boolean atRoot);
1333extern void computeAncestralStates(tree *tr, double referenceLikelihood);
1334extern void makeP_Flex(double z1, double z2, double *rptr, double *EI,  double *EIGN, int numberOfCategories, double *left, double *right, const int numStates);
1335
1336extern void *rax_malloc( size_t size );
1337extern void *rax_realloc(void *p, size_t size, boolean needsMemoryAlignment);
1338extern void rax_free(void *p);
1339extern void *rax_calloc(size_t n, size_t size);
1340
1341#ifdef _WAYNE_MPI
1342
1343extern boolean computeBootStopMPI(tree *tr, char *bootStrapFileName, analdef *adef, double *pearsonAverage);
1344
1345#endif
1346
1347extern void setPartitionMask(tree *tr, int i, boolean *executeModel);
1348extern void resetPartitionMask(tree *tr, boolean *executeModel);
1349#ifdef _USE_PTHREADS
1350
1351extern size_t getContiguousVectorLength(tree *tr);
1352
1353extern void makenewzClassify(tree *tr, int maxiter, double *result, double *z0, double *x1_start, double *x2_start,
1354                             unsigned char *tipX1,  unsigned char *tipX2, int tipCase, boolean *partitionConverged, int insertion);
1355
1356extern void newviewMultiGrain(tree *tr,  double *x1, double *x2, double *x3, int *_ex1, int *_ex2, int *_ex3, unsigned char *_tipX1, unsigned char *_tipX2, 
1357                              int tipCase, double *_pz, double *_qz, int insertion);
1358
1359extern void addTraverseRobIterative(tree *tr, int branchNumber);
1360extern void insertionsParsimonyIterative(tree *tr, int branchNumber);
1361
1362extern void newviewClassify(tree *tr, branchInfo *b, double *z, int insertion);
1363
1364extern double evalCL(tree *tr, double *x2, int *ex2, unsigned char *tip, double *pz, int insertion);
1365
1366extern void testInsertThoroughIterative(tree *tr, int branchNumber);
1367
1368
1369/* parallel MRE stuff */
1370
1371#define MRE_POSSIBLE_CANDIDATE  1
1372#define MRE_EXCLUDED  2
1373#define MRE_ADDED  3   
1374#define SECTION_CONSTANT 1
1375#define MRE_MIN_AMOUNT_JOBS_PER_THREAD 5
1376
1377
1378/* work tags for parallel regions */
1379
1380#define THREAD_NEWVIEW                0
1381#define THREAD_EVALUATE               1
1382#define THREAD_MAKENEWZ               2
1383#define THREAD_MAKENEWZ_FIRST         3
1384#define THREAD_RATE_CATS              4
1385#define THREAD_EVALUATE_VECTOR        7
1386#define THREAD_ALLOC_LIKELIHOOD       8
1387#define THREAD_COPY_RATE_CATS         9
1388#define THREAD_COPY_INVAR             10
1389#define THREAD_COPY_INIT_MODEL        11
1390#define THREAD_FIX_MODEL_INDICES      12
1391#define THREAD_INIT_PARTITION         13
1392#define THREAD_OPT_INVAR              14
1393#define THREAD_OPT_ALPHA              15
1394#define THREAD_OPT_RATE               16
1395#define THREAD_RESET_MODEL            17
1396#define THREAD_COPY_ALPHA             18
1397#define THREAD_COPY_RATES             19
1398#define THREAD_CAT_TO_GAMMA           20
1399#define THREAD_GAMMA_TO_CAT           21
1400#define THREAD_NEWVIEW_MASKED         22
1401#define THREAD_COPY_PARAMS            26
1402#define THREAD_INIT_EPA                     28
1403#define THREAD_GATHER_LIKELIHOOD            29
1404#define THREAD_INSERT_CLASSIFY              30
1405#define THREAD_INSERT_CLASSIFY_THOROUGH     31
1406#define THREAD_GATHER_PARSIMONY             32
1407#define THREAD_PARSIMONY_INSERTIONS         34
1408#define THREAD_PREPARE_EPA_PARSIMONY        35
1409#define THREAD_CLEANUP_EPA_PARSIMONY        36
1410#define THREAD_PREPARE_BIPS_FOR_PRINT       39
1411#define THREAD_MRE_COMPUTE                  40
1412#define THREAD_NEWVIEW_ANCESTRAL            41
1413#define THREAD_GATHER_ANCESTRAL             42
1414#define THREAD_OPT_SCALER                   43
1415#define THREAD_COPY_LG4X_RATES              44
1416#define THREAD_OPT_LG4X_RATES               45
1417
1418
1419/*
1420
1421parallel tree parsing abandoned ...
1422
1423  #define THREAD_FILL_HASH_FOR_CONSENSUS      41
1424
1425parallel dropset comp. currently doesn't scale well
1426
1427  #define THREAD_FIND_BEST_DROPSET            42
1428  #define THREAD_CALC_DROPSETS                43
1429
1430*/
1431
1432/*
1433  Pthreads-based MP computations don't really scale for the
1434  SSE3-based optimized version
1435
1436  #define THREAD_FAST_EVALUATE_PARSIMONY        XX
1437  #define THREAD_FAST_NEWVIEW_PARSIMONY         XX
1438  #define THREAD_INIT_FAST_PARSIMONY            XX
1439*/
1440
1441typedef struct
1442{
1443  tree *tr;
1444  int threadNumber;
1445}
1446  threadData;
1447
1448void threadMakeVector(tree *tr, int tid);
1449void threadComputeAverage(tree *tr, int tid);
1450void threadComputePearson(tree *tr, int tid);
1451extern void optRateCatPthreads(tree *tr, double lower_spacing, double upper_spacing, double *lhs, int n, int tid);
1452extern void masterBarrier(int jobType, tree *tr);
1453
1454#endif
1455
1456boolean isGap(unsigned int *x, int pos);
1457boolean noGap(unsigned int *x, int pos);
1458
1459#ifdef __AVX
1460
1461void newviewGTRGAMMAPROT_AVX_LG4(int tipCase,
1462                                 double *x1, double *x2, double *x3, double *extEV[4], double *tipVector[4],
1463                                 int *ex3, unsigned char *tipX1, unsigned char *tipX2, int n, 
1464                                 double *left, double *right, int *wgt, int *scalerIncrement, const boolean useFastScaling);
1465
1466void newviewGTRCAT_AVX_GAPPED_SAVE(int tipCase,  double *EV,  int *cptr,
1467                                   double *x1_start, double *x2_start,  double *x3_start, double *tipVector,
1468                                   int *ex3, unsigned char *tipX1, unsigned char *tipX2,
1469                                   int n,  double *left, double *right, int *wgt, int *scalerIncrement, const boolean useFastScaling,
1470                                   unsigned int *x1_gap, unsigned int *x2_gap, unsigned int *x3_gap,
1471                                   double *x1_gapColumn, double *x2_gapColumn, double *x3_gapColumn, const int maxCats);
1472
1473void newviewGTRCATPROT_AVX_GAPPED_SAVE(int tipCase, double *extEV,
1474                                       int *cptr,
1475                                       double *x1, double *x2, double *x3, double *tipVector,
1476                                       int *ex3, unsigned char *tipX1, unsigned char *tipX2,
1477                                       int n, double *left, double *right, int *wgt, int *scalerIncrement, const boolean useFastScaling,
1478                                       unsigned int *x1_gap, unsigned int *x2_gap, unsigned int *x3_gap,
1479                                       double *x1_gapColumn, double *x2_gapColumn, double *x3_gapColumn, const int maxCats);
1480
1481void  newviewGTRGAMMA_AVX_GAPPED_SAVE(int tipCase,
1482                                      double *x1_start, double *x2_start, double *x3_start,
1483                                      double *extEV, double *tipVector,
1484                                      int *ex3, unsigned char *tipX1, unsigned char *tipX2,
1485                                      const int n, double *left, double *right, int *wgt, int *scalerIncrement, const boolean useFastScaling,
1486                                      unsigned int *x1_gap, unsigned int *x2_gap, unsigned int *x3_gap, 
1487                                      double *x1_gapColumn, double *x2_gapColumn, double *x3_gapColumn
1488                                      );
1489
1490void newviewGTRGAMMAPROT_AVX_GAPPED_SAVE(int tipCase,
1491                                         double *x1_start, double *x2_start, double *x3_start, double *extEV, double *tipVector,
1492                                         int *ex3, unsigned char *tipX1, unsigned char *tipX2, int n, 
1493                                         double *left, double *right, int *wgt, int *scalerIncrement, const boolean useFastScaling,
1494                                         unsigned int *x1_gap, unsigned int *x2_gap, unsigned int *x3_gap, 
1495                                         double *x1_gapColumn, double *x2_gapColumn, double *x3_gapColumn); 
1496
1497void newviewGTRCAT_AVX(int tipCase,  double *EV,  int *cptr,
1498    double *x1_start, double *x2_start,  double *x3_start, double *tipVector,
1499    int *ex3, unsigned char *tipX1, unsigned char *tipX2,
1500    int n,  double *left, double *right, int *wgt, int *scalerIncrement, const boolean useFastScaling);
1501
1502
1503void newviewGenericCATPROT_AVX(int tipCase, double *extEV,
1504    int *cptr,
1505    double *x1, double *x2, double *x3, double *tipVector,
1506    int *ex3, unsigned char *tipX1, unsigned char *tipX2,
1507    int n, double *left, double *right, int *wgt, int *scalerIncrement, const boolean useFastScaling);
1508
1509
1510void newviewGTRGAMMA_AVX(int tipCase,
1511    double *x1_start, double *x2_start, double *x3_start,
1512    double *EV, double *tipVector,
1513    int *ex3, unsigned char *tipX1, unsigned char *tipX2,
1514    const int n, double *left, double *right, int *wgt, int *scalerIncrement, const boolean useFastScaling);
1515
1516void newviewGTRGAMMAPROT_AVX(int tipCase,
1517                             double *x1, double *x2, double *x3, double *extEV, double *tipVector,
1518                             int *ex3, unsigned char *tipX1, unsigned char *tipX2, int n, 
1519                             double *left, double *right, int *wgt, int *scalerIncrement, const boolean useFastScaling);
1520
1521void newviewGTRCATPROT_AVX(int tipCase, double *extEV,
1522                               int *cptr,
1523                               double *x1, double *x2, double *x3, double *tipVector,
1524                               int *ex3, unsigned char *tipX1, unsigned char *tipX2,
1525                           int n, double *left, double *right, int *wgt, int *scalerIncrement, const boolean useFastScaling);
1526
1527#endif
1528
1529
1530
Note: See TracBrowser for help on using the repository browser.