-
Notifications
You must be signed in to change notification settings - Fork 5
/
bayes.h
1762 lines (1635 loc) · 98.7 KB
/
bayes.h
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
#ifndef __BAYES_H__
#define __BAYES_H__
#include <assert.h>
#include <ctype.h>
#include <float.h>
#include <limits.h>
#include <math.h>
#include <memory.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <stdarg.h>
#include <time.h>
#ifdef OPENMP
#include <omp.h>
#endif
#ifdef USECONFIG_H
# include "config.h"
#elif !defined (XCODE_VERSION) /* some defaults that would otherwise be guessed by configure */
# define PACKAGE_NAME "mrbayes"
# define PACKAGE_VERSION "3.2"
# undef HAVE_LIBREADLINE
# define UNIX_VERSION 1
# define SSE_ENABLED 1
# undef MPI_ENABLED
# undef BEAGLE_ENABLED
# undef FAST_LOG
// chl: add avx related
# undef AVX_ENABLED
# undef FMA_ENABLED
# undef MPI_ENABLED
# undef BEAGLE_ENABLED
#endif
#if defined (MPI_ENABLED)
#include "mpi.h"
#endif
#if defined (BEAGLE_ENABLED)
#include "libhmsbeagle/beagle.h"
#endif
/* Set SSE_ENABLED if SSE SIMD extensions available. */
#ifdef HAVE_SSE
#define SSE_ENABLED
#endif
/* Set AVX_ENABLED if AVX SIMD extensions available. */
#ifdef HAVE_AVX
#define AVX_ENABLED
#endif
/* Set FMA_ENABLED if FMA SIMD extensions available. */
#if defined(HAVE_FMA3) || defined(HAVE_FMA4)
#define FMA_ENABLED
#endif
/* uncomment the following line when releasing, also modify the VERSION_NUMBER below */
#define RELEASE
#ifdef RELEASE
#define VERSION_NUMBER "3.2.6"
#else
#define VERSION_NUMBER "3.2.7-svn"
#endif
#if !defined (UNIX_VERSION) && !defined (WIN_VERSION) && !defined (MAC_VERSION)
# ifdef __MWERKS__
# define MAC_VERSION
# elif defined __APPLE__
# define MAC_VERSION
# else
# define WIN_VERSION
# endif
#endif
#if defined (WIN_VERSION)
# include <windows.h>
# include <winbase.h>
#endif
/* Previous problems with bitfield operations may have been caused by several things. One potential
problem has been that MrBayes has used signed ints or signed longs for the bit operations, and
not all compilers perform bitfield operations on signed ints and longs in the expected manner
for unsigned ints and longs. Actually, the ANSI standard does not specify how to implement
bit operations on signed operands.
Another problem is that some bitshift operations in MrBayes have been performed on an integer
constant "1". This is probably interpreted by most compilers as a signed int and the result of
the bitshift operation is therefore also a signed int. When a signed int has had a different
number of bits than the type specified as SafeLong, this may also have resulted in errors on
some systems.
Both of these problems are fixed now by using unsigned longs for the bitfield operations and
by setting the "1" value needed by some bitfield functions using a variable defined as the same
type as the bitfield containers themselves. Furthermore, I have separated out the cases where
a signed long is required or is more appropriate than the container type of bitfields.
FR 2013-07-06
*/
typedef unsigned long BitsLong;
typedef long RandLong;
#define MRBFLT_MAX DBL_MAX /* maximum possible value that can be stored in MrBFlt */
#define MRBFLT_MIN DBL_MIN /* minimum possible value that can be stored in MrBFlt */
#define MRBFLT_NEG_MAX (-DBL_MAX) /* maximum possible negative value that can be stored in MrBFlt */
typedef double MrBFlt; /* double used for parameter values and generally for floating point values,
if set to float MPI would not work becouse of use MPI_DOUBLE */
typedef float CLFlt; /* single-precision float used for cond likes (CLFlt) to increase speed and reduce memory requirement */
/* set CLFlt to double if you want increased precision */
/* NOTE: CLFlt = double not compatible with SSE_ENABLED */
/*
* Make sure we define SIMD instruction flags in a stepwise manner. That is, if we have FMA, make sure we have AVX;
* if we have AVX, make sure we have SSE.
*/
//#if defined (FMA_ENABLED)
//# if !defined (AVX_ENABLED)
//# define AVX_ENABLED
//# endif
//#endif
//
//#if defined (AVX_ENABLED)
//# if !defined (SSE_ENABLED)
//# define SSE_ENABLED
//# endif
//#endif
/* Define a compiler and vector size for the SSE code */
#if defined (SSE_ENABLED)
//chl: have changed FLOATS_PER_VEC from 4 to 8
# define FLOATS_PER_VEC 8
# if defined (WIN_VERSION)
# define MS_VCPP_SSE
# include <xmmintrin.h>
# else
# define GCC_SSE
# undef ICC_SSE
# include <xmmintrin.h>
# endif
#endif
/* Include header for AVX code */
#if defined (AVX_ENABLED)
# include <immintrin.h>
#endif
#if defined (FMA_ENABLED)
# include <immintrin.h>
#endif
/* Define vector code constants */
//#define VEC_NONE 0
//#define VEC_SSE 1
//#define VEC_AVX 2
//#define VEC_FMA 3
/* For comparing floating points: two values are the same if the absolute difference is less then
this value.
*/
#ifndef ETA
#define ETA (1E-30)
#endif
/* This defined DEBUG() is not used anywhere
#if defined (DEBUGOUTPUT)
#define DEBUG(fmt, arg) printf("%s:%d ",__FILE__,__LINE__); printf(fmt,arg);
#endif
*/
/* TEMPSTRSIZE determines size of temporary sprintf buffer (for SafeSprintf) */
/* A patch was sent in by Allen Smith for SafeSprintf, but I could not get
it compiled on SGI IRIX 6.5 (too old?) with _xpg5_vsnprintf undefined.
The code below is a hack so SafeSprintf never has to reallocate memory.
*/
#ifdef __sgi
#define TEMPSTRSIZE 1000
#else
#define TEMPSTRSIZE 200
#endif
/* NO_ERROR is defined in bayes.h (as 0) and also in WinError.h (as 0L)
ERROR is defined in bayes.h (as 1) and also in WinGDI.h (as 0). we use the bayes.h value */
#undef NO_ERROR
#undef ERROR
#define NO_ERROR 0
#define ERROR 1
#define NO_ERROR_QUIT 2
#define ABORT 3
#define SKIP_COMMAND 4
#undef FALSE
#undef TRUE
#define FALSE 0
#define TRUE 1
#define NO 0
#define YES 1
#define UP 0
#define DOWN 1
#define UPPER 0
#define MIDDLE 1
#define LOWER 2
#define NONINTERACTIVE 0
#define INTERACTIVE 1
#define STANDARD_USER 1
#define DEVELOPER 3
#define DIFFERENT 0
#define SAME 1
#define CONSISTENT_WITH 2
#define LINETERM_UNIX 0
#define LINETERM_MAC 1
#define LINETERM_DOS 2
#define SCREENWIDTH 60
#define SCREENWIDTH2 61
#define AVGSTDDEV 0
#define MAXSTDDEV 1
#define NONE 0
#define DNA 1
#define RNA 2
#define PROTEIN 3
#define RESTRICTION 4
#define STANDARD 5
#define MIXED 6
#define CONTINUOUS 7
#define AAMODEL_POISSON 0
#define AAMODEL_JONES 1
#define AAMODEL_DAY 2
#define AAMODEL_MTREV 3
#define AAMODEL_MTMAM 4
#define AAMODEL_WAG 5
#define AAMODEL_RTREV 6
#define AAMODEL_CPREV 7
#define AAMODEL_VT 8
#define AAMODEL_BLOSUM 9
#define AAMODEL_LG 10
#define AAMODEL_EQ 11
#define AAMODEL_GTR 12 /* aa models with free parameters must be listed last */
#define NUCMODEL_4BY4 0
#define NUCMODEL_DOUBLET 1
#define NUCMODEL_CODON 2
#define NUCMODEL_AA 3
#define NST_MIXED -1 /* anything other than 1, 2, or 6 */
#define MISSING 10000000
#define GAP 10000001
#define UNORD 0
#define ORD 1
#define DOLLO 2
#define IRREV 3
#define IN_CMD 0
#define IN_FILE 1
#define NOTHING 0
#define COMMAND 1
#define PARAMETER 2
#define EQUALSIGN 3
#define COLON 4
#define SEMICOLON 5
#define COMMA 6
#define POUNDSIGN 7
#define QUESTIONMARK 8
#define DASH 9
#define LEFTPAR 10
#define RIGHTPAR 11
#define LEFTCOMMENT 12
#define RIGHTCOMMENT 13
#define ALPHA 14
#define NUMBER 15
#define RETURNSYMBOL 16
#define ASTERISK 17
#define BACKSLASH 18
#define FORWARDSLASH 19
#define EXCLAMATIONMARK 20
#define PERCENT 21
#define QUOTATIONMARK 22
#define WEIRD 23
#define UNKNOWN_TOKEN_TYPE 24
#define LEFTCURL 25
#define RIGHTCURL 26
#define DOLLAR 27
#define AMPERSAND 28
#define VERTICALBAR 29
#define MAX_Q_RATE 100.0f
#define MIN_SHAPE_PARAM 0.00001f
#define MAX_SHAPE_PARAM 100.0f
#define MAX_SITE_RATE 10.0f
#define MAX_GAMMA_CATS 20
#define MAX_GAMMA_CATS_SQUARED 400
#define BRLENS_MIN 0.00000001f // 1E-8f
#define BRLENS_MAX 100.0f
/* BRLENS_MIN must be bigger than TIME_MIN */
#define TIME_MIN 1.0E-11f
#define TIME_MAX 100.0f
#define RELBRLENS_MIN 0.00000001f // 1E-8f
#define RELBRLENS_MAX 100.0f
#define KAPPA_MIN 0.001f
#define KAPPA_MAX 1000.0f
#define GROWTH_MIN 0.000001f
#define GROWTH_MAX 1000000.0f
#define RATE_MIN 0.000001f
#define RATE_MAX 100.0f
#define CPPRATEMULTIPLIER_MIN 0.001f
#define CPPRATEMULTIPLIER_MAX 1000.0f
#define SYMPI_MIN 0.000001f
#define SYMPI_MAX 100.0f
#define ALPHA_MIN 0.0001f
#define ALPHA_MAX 10000.0f
#define DIR_MIN 0.000001f
#define PI_MIN 0.000001f
#define OFFSETEXPLAMBDA_MIN 0.000001f
#define OFFSETEXPLAMBDA_MAX 100000.0f
#define TREEHEIGHT_MIN 0.00000001f
#define TREEHEIGHT_MAX 1000.0f
#define TREEAGE_MIN 0.00000001f
#define TREEAGE_MAX 1000000.0f
#define CPPLAMBDA_MIN 0.00001f
#define CPPLAMBDA_MAX 100.0f
#define TK02VAR_MIN 0.000001f
#define TK02VAR_MAX 10000.0f
#define IGRVAR_MIN 0.000001f
#define IGRVAR_MAX 10000.0f
#define MIXEDVAR_MIN 0.000001f
#define MIXEDVAR_MAX 10000.0f
#define OMEGA_MIN 0.001f
#define OMEGA_MAX 1000.0f
#define POS_MIN 1E-25f
#define POS_MAX 1E25f
#define POS_INFINITY 1E25f
#define NEG_INFINITY -1000000.0f
#define POSREAL_MIN 1E-25f
#define POSREAL_MAX 1E25f
#define CMD_STRING_LENGTH 100000
#define pos(i,j,n) ((i)*(n)+(j))
#define NUM_ALLOCS 100
#define ALLOC_MATRIX 0
#define ALLOC_CHARINFO 2
#define ALLOC_CHARSETS 3
#define ALLOC_TAXA 4
#define ALLOC_TMPSET 5
#define ALLOC_PARTITIONS 6
#define ALLOC_PARTITIONVARS 7
#define ALLOC_TAXASETS 8
#define ALLOC_CONSTRAINTS 9
#define ALLOC_USERTREE 10
#define ALLOC_SUMTPARAMS 11
#define ALLOC_TERMSTATE 12
#define ALLOC_ISPARTAMBIG 13
#define ALLOC_AVAILNODES 25
#define ALLOC_AVAILINDICES 26
#define ALLOC_CURLNL 28
#define ALLOC_CURLNPR 29
#define ALLOC_CHAINID 30
#define ALLOC_PARAMS 31
#define ALLOC_TREE 32
#define ALLOC_NODES 33
#define ALLOC_LOCTAXANAMES 34
#define ALLOC_COMPMATRIX 39
#define ALLOC_NUMSITESOFPAT 40
#define ALLOC_COMPCOLPOS 41
#define ALLOC_COMPCHARPOS 42
#define ALLOC_ORIGCHAR 43
#define ALLOC_PARAMVALUES 46
#define ALLOC_MCMCTREES 47
#define ALLOC_MOVES 48
#define ALLOC_PRELIKES 52
#define ALLOC_SITEJUMP 54
#define ALLOC_MARKOVTIS 55
#define ALLOC_RATEPROBS 56
#define ALLOC_STDTYPE 57
#define ALLOC_PACKEDTREES 58
#define ALLOC_SUMPSTRING 62
#define ALLOC_SUMPINFO 63
#define ALLOC_SWAPINFO 64
#define ALLOC_SYMPIINDEX 65
#define ALLOC_POSSELPROBS 66
#define ALLOC_PBF 68
#define ALLOC_LOCALTAXONCALIBRATION 69
#define ALLOC_SPR_PARSSETS 72
#define ALLOC_PFCOUNTERS 74
#define ALLOC_FILEPOINTERS 75
#define ALLOC_STATS 76
#define ALLOC_DIAGNTREE 77
#define ALLOC_USEDMOVES 82
#define ALLOC_MODEL 83
#define ALLOC_STDSTATEFREQS 84
#define ALLOC_PRINTPARAM 85
#define ALLOC_TREELIST 86
#define ALLOC_TFILEPOS 87
#define ALLOC_BEST 88
#define ALLOC_SPECIESPARTITIONS 89
#define ALLOC_SS 90
#define ALLOC_SAMPLEFOSSILSLICE 91
#define LINKED 0
#define UNLINKED 1
/*paramType*/
#define NUM_LINKED 31
#define P_TRATIO 0
#define P_REVMAT 1
#define P_OMEGA 2
#define P_PI 3
#define P_SHAPE 4
#define P_PINVAR 5
#define P_CORREL 6
#define P_SWITCH 7
#define P_RATEMULT 8
#define P_TOPOLOGY 9
#define P_BRLENS 10
#define P_SPECRATE 11
#define P_EXTRATE 12
#define P_FOSLRATE 13
#define P_POPSIZE 14
#define P_AAMODEL 15
#define P_BRCORR 16
#define P_BRSIGMA 17
#define P_GROWTH 18
#define P_CPPMULTDEV 19
#define P_CPPRATE 20
#define P_CPPEVENTS 21
#define P_TK02VAR 22
#define P_TK02BRANCHRATES 23
#define P_IGRVAR 24
#define P_IGRBRANCHRATES 25
#define P_CLOCKRATE 26
#define P_SPECIESTREE 27
#define P_GENETREERATE 28
#define P_MIXEDVAR 29
#define P_MIXEDBRCHRATES 30
/* NOTE: If you add another parameter, change NUM_LINKED */
// #define CPPm 0 /* CPP rate multipliers */
// #define CPPi 1 /* CPP independent rates */
#define RCL_TK02 0
#define RCL_IGR 1 /* type of mixed relaxed clock model */
#define MAX_NUM_USERTREES 200 /* maximum number of user trees MrBayes will read */
#define MAX_CHAINS 256 /* maximum numbder of chains you can run actually only half of it becouse of m->lnLike[MAX_CHAINS] */
// #define PARAM_NAME_SIZE 400
typedef void * VoidPtr;
typedef int (*CmdFxn)(void);
typedef int (*ParmFxn)(char *, char *);
/* typedef for a ln prior prob fxn */
typedef MrBFlt (*LnPriorProbFxn)(MrBFlt val, MrBFlt *priorParams);
/* typedef for a ln prior prob ratio fxn */
typedef MrBFlt (*LnPriorRatioFxn)(MrBFlt newVal, MrBFlt oldVal, MrBFlt *priorParams);
typedef struct
{
MrBFlt sum; /* sum of standard deviations */
MrBFlt max; /* maximum standard deviation */
MrBFlt numPartitions;
MrBFlt numSamples;
MrBFlt avgStdDev;
MrBFlt **pair;
} STATS;
/* enumeration for calibration prior type */
enum CALPRIOR
{
unconstrained,
fixed,
uniform,
offsetExponential,
truncatedNormal,
logNormal,
offsetLogNormal,
standardGamma,
offsetGamma
};
enum ConstraintType
{
PARTIAL,
NEGATIVE,
HARD
};
enum CodingType
{
ALL = 0,
NOABSENCESITES = 1,
NOPRESENCESITES = 2,
VARIABLE = 3,
NOSINGLETONPRESENCE = 4,
NOSINGLETONABSENCE = 8,
NOSINGLETONS = 12,
INFORMATIVE = 15
};
/* typedef for calibration */
typedef struct calibration
{
char name[100];
enum CALPRIOR prior;
MrBFlt priorParams[3];
LnPriorProbFxn LnPriorProb;
LnPriorRatioFxn LnPriorRatio;
MrBFlt min;
MrBFlt max;
}
Calibration;
/* typedef for tree (topology) list element */
typedef struct element
{
struct element *next;
int *order;
} TreeListElement;
/* typedef for list of trees (topologies) */
typedef struct
{
TreeListElement *first;
TreeListElement *last;
} TreeList;
/* typedef for packed tree */
typedef struct
{
int *order;
MrBFlt *brlens;
} PackedTree;
/* typedef for binary tree node */
/* NOTE: Any variable added here must also be copied in CopyTrees */
typedef struct node
{
char *label; /*!< name of node if tip */
struct node *left, *right, *anc; /*!< pointers to adjacent nodes */
int memoryIndex; /*!< memory index (do not change) */
int index; /*!< index to node (0 to numLocalTaxa for tips) */
int upDateCl; /*!< cond likes need update? */
int upDateTi; /*!< transition probs need update? */
int scalerNode; /*!< is node scaling cond likes? */
int isLocked; /*!< is node locked? */
int lockID; /*!< id of lock */
int isDated; /*!< is node dated (calibrated)? */
int marked, x, y; /*!< scratch variables */
MrBFlt d; /*!< scratch variable */
BitsLong *partition; /*!< pointer to bitfield describing splits */
MrBFlt length; /*!< length of pending branch */
MrBFlt nodeDepth; /*!< node depth (height) */
MrBFlt age; /*!< age of node */
Calibration *calibration; /*!< pointer to calibration data */
}
TreeNode;
/* typedef for binary tree */
typedef struct
{
char name[100]; /*!< name of tree */
int memNodes; /*!< number of allocated nodes (do not exceed!) */
int nNodes; /*!< number of nodes in tree (including lower root in rooted trees) */
int nIntNodes; /*!< number of interior nodes in tree (excluding lower root in rooted trees) */
int isRooted; /*!< is tree rooted? */
int isClock; /*!< is tree clock? */
int isCalibrated; /*!< is tree calibrated? */
int nRelParts; /*!< number of relevant partitions */
int *relParts; /*!< pointer to relevant partitions */
int checkConstraints; /*!< does tree have constraints? */
int nConstraints; /*!< number of constraints */
int *constraints; /*!< pointer to constraints */
int nLocks; /*!< number of constrained (locked) nodes */
TreeNode **allDownPass; /*!< downpass array of all nodes */
TreeNode **intDownPass; /*!< downpass array of interior nodes (including upper but excluding lower root in rooted trees) */
TreeNode *root; /*!< pointer to root (lower root in rooted trees) */
TreeNode *nodes; /*!< array containing the nodes */
BitsLong *bitsets; /*!< pointer to bitsets describing splits */
BitsLong *flags; /*!< pointer to cond like flags */
int fromUserTree; /*!< YES is set for the trees whoes branch lengthes are set from user tree(as start tree or fix branch length prior), NO otherwise */
}
Tree;
/* typedef for node in polytomous tree */
typedef struct pNode
{
char label[100]; /*!< name of node if terminal */
struct pNode *left, *sib, *anc; /*!< pointers to adjacent nodes */
int x, y, mark; /*!< scratch variables */
int partitionIndex; /*!< partition index in sumt (scratch) */
int index; /*!< index of node (if < numLocalTaxa = local taxon index) */
int memoryIndex; /*!< immutable index of memory position */
int isLocked; /*!< is the node locked? */
int lockID; /*!< id of lock */
int isDated; /*!< is node dated? */
MrBFlt length; /*!< age of node */
MrBFlt depth; /*!< depth (height) of node */
MrBFlt age; /*!< age of node */
MrBFlt support, f; /*!< scratch variables */
BitsLong *partition; /*!< pointer to partition (split) bitset */
Calibration *calibration; /*!< pointer to dating of node */
}
PolyNode;
/* typedef for polytomous tree */
typedef struct
{
char name[100]; /*!< name of tree */
int memNodes; /*!< number of allocated nodes; do not exceed! */
int nNodes; /*!< number of nodes in tree */
int nIntNodes; /*!< number of interior nodes in tree */
PolyNode **allDownPass; /*!< downpass array over all nodes */
PolyNode **intDownPass; /*!< downpass array over interior nodes */
PolyNode *root; /*!< pointer to root (lower for rooted trees */
PolyNode *nodes; /*!< array holding the tree nodes */
BitsLong *bitsets; /*!< bits describing partitions (splits) */
int nBSets; /*!< number of effective branch length sets */
int nESets; /*!< number of breakpoint rate sets */
char **bSetName; /*!< names of effective branch length sets */
char **eSetName; /*!< names of breakpoint rate sets */
int **nEvents; /*!< number of branch events of bp rate set */
MrBFlt ***position; /*!< position of branch events */
MrBFlt ***rateMult; /*!< parameter of branch events */
MrBFlt **effectiveBrLen; /*!< effective branch lengths of ebl set */
int brlensDef; /*!< are brlens defined ? */
int isRooted; /*!< is tree rooted? */
int isClock; /*!< is tree clock? */
int isCalibrated; /*!< is tree calibrated? */
int isRelaxed; /*!< is tree relaxed? */
MrBFlt clockRate; /*!< clock rate */
int popSizeSet; /*!< does tree have a population size set? */
MrBFlt *popSize; /*!< the population size */
char *popSizeSetName; /*!< name of the population size set */
}
PolyTree;
/* struct for holding model parameter info for the mcmc run */
typedef struct param
{
int index; /* index to the parameter (0, 1, 2, ...) */
int paramType; /* the type of the parameter */
int paramId; /* unique ID for parameter x prior combination */
MrBFlt *values; /* main values of parameter */
MrBFlt *subValues; /* subvalues of parameter */
int *intValues; /* integer values (model index/growth fxn) */
int nValues; /* number of values */
int nSubValues; /* number of subvalues */
int nIntValues; /* number of intvalues */
MrBFlt min; /* minimum value of parameter */
MrBFlt max; /* maximum value of parameter */
int *relParts; /* pointer to relevant divisions */
int nRelParts; /* number of relevant divisions */
int upDate; /* update flag (for copying) */
struct param **subParams; /* pointers to subparams (for topology) */
int nSubParams; /* number of subparams */
Tree **tree; /* pointer to tree ptrs (for brlens & topology) */
int treeIndex; /* index to first tree in mcmcTree */
int hasBinaryStd; /* has binary standard chars */
int *sympiBsIndex; /* pointer to sympi bsIndex (std chars) */
int *sympinStates; /* pointer to sympi nStates (std chars) */
int *sympiCType; /* pointer to sympi cType (std chars) */
int nSympi; /* number of sympis */
int printParam; /* whether parameter should be printed */
int nPrintSubParams; /* number of subparams that should be printed */
char *paramHeader; /* a string holding header for param values */
char *name; /* string holding name of parameter */
char *paramTypeName; /* pointer to description of parameter type */
int checkConstraints; /* is tree parameter constrained? */
int fill; /* flags whether the parameter should be filled */
int nStdStateFreqs; /* number of std state frequencies */
MrBFlt *stdStateFreqs; /* pointer to std state frequencies */
int **nEvents; /* number of branch events for Cpp model */
/* nEvents[0..2*numCains][0..numNodes=2*numTaxa] */
MrBFlt ***position; /* event positions for Cpp relaxed clock model */
MrBFlt ***rateMult; /* rate multipliers for Cpp relaxed clock model */
int affectsLikelihood; /* does parameter directly influence likelihood? */
MrBFlt* priorParams; /* pointer to the prior parameters */
LnPriorProbFxn LnPriorProb; /* ln prior prob function */
LnPriorRatioFxn LnPriorRatio; /* ln prior prob ratio function */
} Param;
#if defined(THREADS_ENABLED)
#include <pthread.h>
typedef struct s_launch_struct
{
int chain;
int division;
MrBFlt* lnL;
} LaunchStruct;
#endif
/* parameter ID values */
/* identifies unique model parameter x prior combinations */
#define TRATIO_DIR 1
#define TRATIO_FIX 2
#define REVMAT_DIR 3
#define REVMAT_FIX 4
#define OMEGA_DIR 5
#define OMEGA_FIX 6
#define SYMPI_UNI 7
#define SYMPI_UNI_MS 8
#define SYMPI_EXP 9
#define SYMPI_EXP_MS 10
#define SYMPI_FIX 11
#define SYMPI_FIX_MS 12
#define SYMPI_EQUAL 13
#define PI_DIR 14
#define PI_USER 15
#define PI_EMPIRICAL 16
#define PI_EQUAL 17
#define PI_FIXED 18
#define SHAPE_UNI 19
#define SHAPE_EXP 20
#define SHAPE_FIX 21
#define PINVAR_UNI 22
#define PINVAR_FIX 23
#define CORREL_UNI 24
#define CORREL_FIX 25
#define SWITCH_UNI 26
#define SWITCH_EXP 27
#define SWITCH_FIX 28
#define RATEMULT_DIR 29
#define RATEMULT_FIX 30
#define TOPOLOGY_NCL_UNIFORM 31
#define TOPOLOGY_NCL_CONSTRAINED 32
#define TOPOLOGY_NCL_FIXED 33
#define TOPOLOGY_NCL_UNIFORM_HOMO 34
#define TOPOLOGY_NCL_UNIFORM_HETERO 35
#define TOPOLOGY_NCL_CONSTRAINED_HOMO 36
#define TOPOLOGY_NCL_CONSTRAINED_HETERO 37
#define TOPOLOGY_NCL_FIXED_HOMO 38
#define TOPOLOGY_NCL_FIXED_HETERO 39
#define TOPOLOGY_CL_UNIFORM 40
#define TOPOLOGY_CL_CONSTRAINED 41
#define TOPOLOGY_CL_FIXED 42
#define TOPOLOGY_CCL_UNIFORM 43
#define TOPOLOGY_CCL_CONSTRAINED 44
#define TOPOLOGY_CCL_FIXED 45
#define TOPOLOGY_PARSIMONY_UNIFORM 46
#define TOPOLOGY_PARSIMONY_CONSTRAINED 47
#define TOPOLOGY_PARSIMONY_FIXED 48
#define BRLENS_UNI 49
#define BRLENS_EXP 50
#define BRLENS_GamDir 51
#define BRLENS_iGmDir 52
#define BRLENS_twoExp 53
#define BRLENS_FIXED 54
#define BRLENS_CLOCK_UNI 55
#define BRLENS_CLOCK_COAL 56
#define BRLENS_CLOCK_BD 57
#define BRLENS_CLOCK_FIXED 58
#define BRLENS_CLOCK_SPCOAL 59
#define BRLENS_CLOCK_FOSSIL 60
#define BRLENS_PARSIMONY 61
#define SPECRATE_UNI 62
#define SPECRATE_EXP 63
#define SPECRATE_FIX 64
#define EXTRATE_BETA 65
#define EXTRATE_FIX 66
#define FOSLRATE_BETA 67
#define FOSLRATE_FIX 68
#define POPSIZE_UNI 69
#define POPSIZE_GAMMA 70
#define POPSIZE_FIX 71
#define POPSIZE_NORMAL 72
#define POPSIZE_LOGNORMAL 73
#define AAMODEL_FIX 74
#define AAMODEL_MIX 75
#define GROWTH_UNI 76
#define GROWTH_EXP 77
#define GROWTH_FIX 78
#define GROWTH_NORMAL 79
#define OMEGA_BUD 80
#define OMEGA_BUF 81
#define OMEGA_BED 82
#define OMEGA_BEF 83
#define OMEGA_BFD 84
#define OMEGA_BFF 85
#define OMEGA_FUD 86
#define OMEGA_FUF 87
#define OMEGA_FED 88
#define OMEGA_FEF 89
#define OMEGA_FFD 90
#define OMEGA_FFF 91
#define OMEGA_ED 92
#define OMEGA_EF 93
#define OMEGA_FD 94
#define OMEGA_FF 95
#define OMEGA_10UUB 96
#define OMEGA_10UUF 97
#define OMEGA_10UEB 98
#define OMEGA_10UEF 99
#define OMEGA_10UFB 100
#define OMEGA_10UFF 101
#define OMEGA_10EUB 102
#define OMEGA_10EUF 103
#define OMEGA_10EEB 104
#define OMEGA_10EEF 105
#define OMEGA_10EFB 106
#define OMEGA_10EFF 107
#define OMEGA_10FUB 108
#define OMEGA_10FUF 109
#define OMEGA_10FEB 110
#define OMEGA_10FEF 111
#define OMEGA_10FFB 112
#define OMEGA_10FFF 113
#define CPPRATE_FIX 114
#define CPPRATE_EXP 115
#define CPPMULTDEV_FIX 116
#define TK02VAR_FIX 117
#define TK02VAR_EXP 118
#define TK02VAR_UNI 119
#define TOPOLOGY_RCL_UNIFORM 120
#define TOPOLOGY_RCL_CONSTRAINED 121
#define TOPOLOGY_RCL_FIXED 122
#define TOPOLOGY_RCCL_UNIFORM 123
#define TOPOLOGY_RCCL_CONSTRAINED 124
#define TOPOLOGY_RCCL_FIXED 125
#define TOPOLOGY_SPECIESTREE 126
#define CPPEVENTS 127
#define TK02BRANCHRATES 128
#define TOPOLOGY_FIXED 129
#define IGRVAR_FIX 130
#define IGRVAR_EXP 131
#define IGRVAR_UNI 132
#define IGRBRANCHRATES 133
#define CLOCKRATE_FIX 134
#define CLOCKRATE_NORMAL 135
#define CLOCKRATE_LOGNORMAL 136
#define CLOCKRATE_GAMMA 137
#define CLOCKRATE_EXP 138
#define SPECIESTREE_UNIFORM 139
#define GENETREERATEMULT_DIR 140
#define GENETREERATEMULT_FIX 141
#define REVMAT_MIX 142
#define MIXEDVAR_FIX 143
#define MIXEDVAR_EXP 144
#define MIXEDVAR_UNI 145
#define MIXEDBRCHRATES 146
#if defined (BEAGLE_ENABLED)
#define MB_BEAGLE_SCALE_ALWAYS 0
#define MB_BEAGLE_SCALE_DYNAMIC 1
#if defined (_DEBUG)
#define MB_PRINT_DYNAMIC_RESCALE_FAIL_STAT
#endif
#endif
/* typedef for a MoveFxn */
typedef int (MoveFxn)(Param *param, int chain, RandLong *seed, MrBFlt *lnPriorRatio, MrBFlt *lnProposalRatio, MrBFlt *mvp);
/* typedef for an ApplicFxn */
typedef int (ApplicFxn)(Param *param);
/* typedef for an AutotuneFxn */
typedef void (AutotuneFxn)(MrBFlt acceptanceRate, MrBFlt targetRate, int batch, MrBFlt *tuningParameter, MrBFlt minTuning, MrBFlt maxTuning);
/* struct holding info on each move type that the program handles */
typedef struct
{
MoveFxn *moveFxn; /* pointer to the move function */
ApplicFxn *isApplicable; /* pointer to function determining whether move is applicable to a parameter */
int nApplicable; /* number of relevant params */
int applicableTo[40]; /* pointer to ID of relevant params */
char *name; /* name of the move type */
char *shortName; /* abbreviated name of the move type */
char *paramName; /* name of subparameter if complex parameter */
int subParams; /* are we changing subparams (brlens of topol.) */
char *tuningName[5]; /* name of tuning params */
char *shortTuningName[5];/* short name of tuning params */
MrBFlt relProposalProb; /* default relative proposal probability */
int numTuningParams; /* number of tuning parameters */
MrBFlt tuningParam[5]; /* default tuning parameters for the proposal */
MrBFlt minimum[5]; /* minimum values for tuning params */
MrBFlt maximum[5]; /* maximum values for tuning params */
int parsimonyBased; /* this move is based on parsimony (YES/NO) */
int level; /* user level of this move */
AutotuneFxn *Autotune; /* pointer to the autotuning function */
MrBFlt targetRate; /* default target acceptance rate for autotuning */
} MoveType;
/* max number of move types */
#define NUM_MOVE_TYPES 100
/* struct holding info on each move */
/* Note: This allows several different moves to affect the same parameter */
/* It also allows the same move to affect different parameters as before */
/* This is also a good place to keep the proposal probs */
typedef struct
{
char *name; /* name of the move */
MoveType *moveType; /* pointer to the move type */
MoveFxn *moveFxn; /* pointer to the move function */
Param *parm; /* ptr to parameter the move applies to */
MrBFlt *relProposalProb; /* the relative proposal probability */
MrBFlt *cumProposalProb; /* the cumulative proposal probability */
int *nAccepted; /* number of accepted moves */
int *nTried; /* number of tried moves */
int *nBatches; /* counter for autotuning rounds */
int *nTotAccepted; /* total number of accepted moves */
int *nTotTried; /* total number of tried moves */
MrBFlt *targetRate; /* target acceptance rate for autotuning */
MrBFlt *lastAcceptanceRate;/* acceptance rate in last complete batch */
MrBFlt **tuningParam; /* tuning parameters for the move */
} MCMCMove;
typedef int (*LikeDownFxn)(TreeNode *, int, int);
typedef int (*LikeRootFxn)(TreeNode *, int, int);
typedef int (*LikeScalerFxn)(TreeNode *, int, int);
typedef int (*LikeFxn)(TreeNode *, int, int, MrBFlt *, int);
typedef int (*TiProbFxn)(TreeNode *, int, int);
typedef int (*LikeUpFxn)(TreeNode *, int, int);
typedef int (*PrintAncStFxn)(TreeNode *, int, int);
typedef int (*StateCodeFxn) (int);
typedef int (*PrintSiteRateFxn) (TreeNode *, int, int);
typedef int (*PosSelProbsFxn) (TreeNode *, int, int);
typedef int (*SiteOmegasFxn) (TreeNode *, int, int);
typedef struct cmdtyp
{
int cmdNumber;
char *string;
int specialCmd;
CmdFxn cmdFxnPtr;
short numParms;
short parmList[50];
int expect;
char *cmdDescription;
int cmdUse;
int hiding;
}
CmdType;
typedef struct parm
{
char *string; /* parameter name */
char *valueList; /* list of values that could be input */
ParmFxn fp; /* function pointer */
}
ParmInfo, *ParmInfoPtr;
typedef struct model
{
int dataType; /* data type for partition */
int nStates; /* number of states for this type of data */
int codon[64]; /* gives protein ID for each codon */
int codonNucs[64][3]; /* gives triplet for each codon */
int codonAAs[64]; /* gives protein ID for implemented code */
char nucModel[100]; /* nucleotide model used */
char nst[100]; /* number of substitution types */
char parsModel[100]; /* use the (so-called) parsimony model */
char geneticCode[100]; /* genetic code used */
int coding; /* type of patterns encoded */
char codingString[100]; /* string describing type of patterns encoded */
char ploidy[100]; /* ploidy level */
char omegaVar[100]; /* type of omega variation model */
char ratesModel[100]; /* rates across sites model */
int numGammaCats; /* number of categories for gamma approximation */
char useGibbs[100]; /* flags whether Gibbs sampling of discrete gamma is used */
int gibbsFreq; /* frequency of Gibbs resampling of discrete gamma */
int numBetaCats; /* number of categories for beta approximation */
int numM10GammaCats; /* number of cats for gamma approx (M10 model) */
int numM10BetaCats; /* number of cats for beta approx (M10 model) */
char covarionModel[100];/* use covarion model? (yes/no) */
char augmentData[100]; /* should data be augmented */
char tRatioPr[100]; /* prior for ti/tv rate ratio */
MrBFlt tRatioFix;
MrBFlt tRatioDir[2];
char revMatPr[100]; /* prior for GTR model */
MrBFlt revMatFix[6];
MrBFlt revMatDir[6];
MrBFlt revMatSymDir; /* prior for mixed GTR subspace model */
char aaModelPr[100]; /* prior for amino acid model */
char aaModel[100];
MrBFlt aaModelPrProbs[10];
char aaRevMatPr[100]; /* prior for aa GTR model */
MrBFlt aaRevMatFix[190];