18 #elif !defined (XCODE_VERSION) /* some defaults that would otherwise be guessed by configure */
19 # define PACKAGE_NAME "mrbayes"
20 # define PACKAGE_VERSION "3.2"
21 # undef HAVE_LIBREADLINE
22 # define UNIX_VERSION 1
23 # define SSE_ENABLED 1
25 # undef BEAGLE_ENABLED
29 #if defined (MPI_ENABLED)
33 #if defined (BEAGLE_ENABLED)
34 #include "libhmsbeagle/beagle.h"
37 /* uncomment the following line when releasing, also modify the VERSION_NUMBER below */
40 #define VERSION_NUMBER "3.2.6"
42 #define VERSION_NUMBER "3.2.7-svn"
45 #if !defined (UNIX_VERSION) && !defined (WIN_VERSION) && !defined (MAC_VERSION)
48 # elif defined __APPLE__
55 #if defined (WIN_VERSION)
60 /* Previous problems with bitfield operations may have been caused by several things. One potential
61 problem has been that MrBayes has used signed ints or signed longs for the bit operations, and
62 not all compilers perform bitfield operations on signed ints and longs in the expected manner
63 for unsigned ints and longs. Actually, the ANSI standard does not specify how to implement
64 bit operations on signed operands.
66 Another problem is that some bitshift operations in MrBayes have been performed on an integer
67 constant "1". This is probably interpreted by most compilers as a signed int and the result of
68 the bitshift operation is therefore also a signed int. When a signed int has had a different
69 number of bits than the type specified as SafeLong, this may also have resulted in errors on
72 Both of these problems are fixed now by using unsigned longs for the bitfield operations and
73 by setting the "1" value needed by some bitfield functions using a variable defined as the same
74 type as the bitfield containers themselves. Furthermore, I have separated out the cases where
75 a signed long is required or is more appropriate than the container type of bitfields.
79 typedef unsigned long BitsLong;
80 typedef long RandLong;
82 #define MRBFLT_MAX DBL_MAX /* maximum possible value that can be stored in MrBFlt */
83 #define MRBFLT_MIN DBL_MIN /* minimum possible value that can be stored in MrBFlt */
84 #define MRBFLT_NEG_MAX (-DBL_MAX) /* maximum possible negative value that can be stored in MrBFlt */
85 typedef double MrBFlt; /* double used for parameter values and generally for floating point values,
86 if set to float MPI would not work becouse of use MPI_DOUBLE */
87 typedef float CLFlt; /* single-precision float used for cond likes (CLFlt) to increase speed and reduce memory requirement */
88 /* set CLFlt to double if you want increased precision */
89 /* NOTE: CLFlt = double not compatible with SSE_ENABLED */
91 /* Define a compiler and vector size for the SSE code */
92 #if defined (SSE_ENABLED)
93 # define FLOATS_PER_VEC 4
94 # if defined (WIN_VERSION)
96 # include <xmmintrin.h>
100 # include <xmmintrin.h>
104 /* For comparing floating points: two values are the same if the absolute difference is less then
111 /* This defined DEBUG() is not used anywhere
112 #if defined (DEBUGOUTPUT)
113 #define DEBUG(fmt, arg) printf("%s:%d ",__FILE__,__LINE__); printf(fmt,arg);
117 /* TEMPSTRSIZE determines size of temporary sprintf buffer (for SafeSprintf) */
118 /* A patch was sent in by Allen Smith for SafeSprintf, but I could not get
119 it compiled on SGI IRIX 6.5 (too old?) with _xpg5_vsnprintf undefined.
120 The code below is a hack so SafeSprintf never has to reallocate memory.
123 #define TEMPSTRSIZE 1000
125 #define TEMPSTRSIZE 200
128 /* NO_ERROR is defined in bayes.h (as 0) and also in WinError.h (as 0L)
129 ERROR is defined in bayes.h (as 1) and also in WinGDI.h (as 0). we use the bayes.h value */
134 #define NO_ERROR_QUIT 2
136 #define SKIP_COMMAND 4
153 #define NONINTERACTIVE 0
154 #define INTERACTIVE 1
156 #define STANDARD_USER 1
161 #define CONSISTENT_WITH 2
163 #define LINETERM_UNIX 0
164 #define LINETERM_MAC 1
165 #define LINETERM_DOS 2
167 #define SCREENWIDTH 60
168 #define SCREENWIDTH2 61
177 #define RESTRICTION 4
182 #define AAMODEL_POISSON 0
183 #define AAMODEL_JONES 1
184 #define AAMODEL_DAY 2
185 #define AAMODEL_MTREV 3
186 #define AAMODEL_MTMAM 4
187 #define AAMODEL_WAG 5
188 #define AAMODEL_RTREV 6
189 #define AAMODEL_CPREV 7
191 #define AAMODEL_BLOSUM 9
192 #define AAMODEL_LG 10
193 #define AAMODEL_EQ 11
194 #define AAMODEL_GTR 12 /* aa models with free parameters must be listed last */
196 #define NUCMODEL_4BY4 0
197 #define NUCMODEL_DOUBLET 1
198 #define NUCMODEL_CODON 2
199 #define NUCMODEL_AA 3
201 #define NST_MIXED -1 /* anything other than 1, 2, or 6 */
203 #define MISSING 10000000
222 #define QUESTIONMARK 8
226 #define LEFTCOMMENT 12
227 #define RIGHTCOMMENT 13
230 #define RETURNSYMBOL 16
233 #define FORWARDSLASH 19
234 #define EXCLAMATIONMARK 20
236 #define QUOTATIONMARK 22
238 #define UNKNOWN_TOKEN_TYPE 24
243 #define VERTICALBAR 29
245 #define MAX_Q_RATE 100.0f
246 #define MIN_SHAPE_PARAM 0.00001f
247 #define MAX_SHAPE_PARAM 100.0f
248 #define MAX_SITE_RATE 10.0f
249 #define MAX_GAMMA_CATS 20
250 #define MAX_GAMMA_CATS_SQUARED 400
251 #define BRLENS_MIN 0.00000001f // 1E-8f
252 #define BRLENS_MAX 100.0f
253 /* BRLENS_MIN must be bigger than TIME_MIN */
254 #define TIME_MIN 1.0E-11f
255 #define TIME_MAX 100.0f
256 #define RELBRLENS_MIN 0.00000001f // 1E-8f
257 #define RELBRLENS_MAX 100.0f
258 #define KAPPA_MIN 0.001f
259 #define KAPPA_MAX 1000.0f
260 #define GROWTH_MIN 0.000001f
261 #define GROWTH_MAX 1000000.0f
262 #define RATE_MIN 0.000001f
263 #define RATE_MAX 100.0f
264 #define CPPRATEMULTIPLIER_MIN 0.001f
265 #define CPPRATEMULTIPLIER_MAX 1000.0f
266 #define SYMPI_MIN 0.000001f
267 #define SYMPI_MAX 100.0f
268 #define ALPHA_MIN 0.0001f
269 #define ALPHA_MAX 10000.0f
270 #define DIR_MIN 0.000001f
271 #define PI_MIN 0.000001f
272 #define OFFSETEXPLAMBDA_MIN 0.000001f
273 #define OFFSETEXPLAMBDA_MAX 100000.0f
274 #define TREEHEIGHT_MIN 0.00000001f
275 #define TREEHEIGHT_MAX 1000.0f
276 #define TREEAGE_MIN 0.00000001f
277 #define TREEAGE_MAX 1000000.0f
278 #define CPPLAMBDA_MIN 0.00001f
279 #define CPPLAMBDA_MAX 100.0f
280 #define TK02VAR_MIN 0.000001f
281 #define TK02VAR_MAX 10000.0f
282 #define IGRVAR_MIN 0.000001f
283 #define IGRVAR_MAX 10000.0f
284 #define MIXEDVAR_MIN 0.000001f
285 #define MIXEDVAR_MAX 10000.0f
286 #define OMEGA_MIN 0.001f
287 #define OMEGA_MAX 1000.0f
289 #define POS_MIN 1E-25f
290 #define POS_MAX 1E25f
291 #define POS_INFINITY 1E25f
292 #define NEG_INFINITY -1000000.0f
293 #define POSREAL_MIN 1E-25f
294 #define POSREAL_MAX 1E25f
296 #define CMD_STRING_LENGTH 100000
298 #define pos(i,j,n) ((i)*(n)+(j))
300 #define NUM_ALLOCS 100
302 #define ALLOC_MATRIX 0
303 #define ALLOC_CHARINFO 2
304 #define ALLOC_CHARSETS 3
306 #define ALLOC_TMPSET 5
307 #define ALLOC_PARTITIONS 6
308 #define ALLOC_PARTITIONVARS 7
309 #define ALLOC_TAXASETS 8
310 #define ALLOC_CONSTRAINTS 9
311 #define ALLOC_USERTREE 10
312 #define ALLOC_SUMTPARAMS 11
313 #define ALLOC_TERMSTATE 12
314 #define ALLOC_ISPARTAMBIG 13
315 #define ALLOC_AVAILNODES 25
316 #define ALLOC_AVAILINDICES 26
317 #define ALLOC_CURLNL 28
318 #define ALLOC_CURLNPR 29
319 #define ALLOC_CHAINID 30
320 #define ALLOC_PARAMS 31
321 #define ALLOC_TREE 32
322 #define ALLOC_NODES 33
323 #define ALLOC_LOCTAXANAMES 34
324 #define ALLOC_COMPMATRIX 39
325 #define ALLOC_NUMSITESOFPAT 40
326 #define ALLOC_COMPCOLPOS 41
327 #define ALLOC_COMPCHARPOS 42
328 #define ALLOC_ORIGCHAR 43
329 #define ALLOC_PARAMVALUES 46
330 #define ALLOC_MCMCTREES 47
331 #define ALLOC_MOVES 48
332 #define ALLOC_PRELIKES 52
333 #define ALLOC_SITEJUMP 54
334 #define ALLOC_MARKOVTIS 55
335 #define ALLOC_RATEPROBS 56
336 #define ALLOC_STDTYPE 57
337 #define ALLOC_PACKEDTREES 58
338 #define ALLOC_SUMPSTRING 62
339 #define ALLOC_SUMPINFO 63
340 #define ALLOC_SWAPINFO 64
341 #define ALLOC_SYMPIINDEX 65
342 #define ALLOC_POSSELPROBS 66
344 #define ALLOC_LOCALTAXONCALIBRATION 69
345 #define ALLOC_SPR_PARSSETS 72
346 #define ALLOC_PFCOUNTERS 74
347 #define ALLOC_FILEPOINTERS 75
348 #define ALLOC_STATS 76
349 #define ALLOC_DIAGNTREE 77
350 #define ALLOC_USEDMOVES 82
351 #define ALLOC_MODEL 83
352 #define ALLOC_STDSTATEFREQS 84
353 #define ALLOC_PRINTPARAM 85
354 #define ALLOC_TREELIST 86
355 #define ALLOC_TFILEPOS 87
356 #define ALLOC_BEST 88
357 #define ALLOC_SPECIESPARTITIONS 89
359 #define ALLOC_SAMPLEFOSSILSLICE 91
365 #define NUM_LINKED 31
377 #define P_SPECRATE 11
379 #define P_FOSLRATE 13
385 #define P_CPPMULTDEV 19
387 #define P_CPPEVENTS 21
389 #define P_TK02BRANCHRATES 23
391 #define P_IGRBRANCHRATES 25
392 #define P_CLOCKRATE 26
393 #define P_SPECIESTREE 27
394 #define P_GENETREERATE 28
395 #define P_MIXEDVAR 29
396 #define P_MIXEDBRCHRATES 30
397 /* NOTE: If you add another parameter, change NUM_LINKED */
399 // #define CPPm 0 /* CPP rate multipliers */
400 // #define CPPi 1 /* CPP independent rates */
402 #define RCL_IGR 1 /* type of mixed relaxed clock model */
404 #define MAX_NUM_USERTREES 200 /* maximum number of user trees MrBayes will read */
405 #define MAX_CHAINS 256 /* maximum numbder of chains you can run actually only half of it becouse of m->lnLike[MAX_CHAINS] */
407 // #define PARAM_NAME_SIZE 400
409 typedef void * VoidPtr;
410 typedef int (*CmdFxn)(void);
411 typedef int (*ParmFxn)(char *, char *);
413 /* typedef for a ln prior prob fxn */
414 typedef MrBFlt (*LnPriorProbFxn)(MrBFlt val, MrBFlt *priorParams);
416 /* typedef for a ln prior prob ratio fxn */
417 typedef MrBFlt (*LnPriorRatioFxn)(MrBFlt newVal, MrBFlt oldVal, MrBFlt *priorParams);
421 MrBFlt sum; /* sum of standard deviations */
422 MrBFlt max; /* maximum standard deviation */
423 MrBFlt numPartitions;
429 /* enumeration for calibration prior type */
456 NOSINGLETONPRESENCE = 4,
457 NOSINGLETONABSENCE = 8,
462 /* typedef for calibration */
463 typedef struct calibration
467 MrBFlt priorParams[3];
468 LnPriorProbFxn LnPriorProb;
469 LnPriorRatioFxn LnPriorRatio;
475 /* typedef for tree (topology) list element */
476 typedef struct element
478 struct element *next;
482 /* typedef for list of trees (topologies) */
485 TreeListElement *first;
486 TreeListElement *last;
489 /* typedef for packed tree */
496 /* typedef for binary tree node */
497 /* NOTE: Any variable added here must also be copied in CopyTrees */
500 char *label; /*!< name of node if tip */
501 struct node *left, *right, *anc; /*!< pointers to adjacent nodes */
502 int memoryIndex; /*!< memory index (do not change) */
503 int index; /*!< index to node (0 to numLocalTaxa for tips) */
504 int upDateCl; /*!< cond likes need update? */
505 int upDateTi; /*!< transition probs need update? */
506 int scalerNode; /*!< is node scaling cond likes? */
507 int isLocked; /*!< is node locked? */
508 int lockID; /*!< id of lock */
509 int isDated; /*!< is node dated (calibrated)? */
510 int marked, x, y; /*!< scratch variables */
511 MrBFlt d; /*!< scratch variable */
512 BitsLong *partition; /*!< pointer to bitfield describing splits */
513 MrBFlt length; /*!< length of pending branch */
514 MrBFlt nodeDepth; /*!< node depth (height) */
515 MrBFlt age; /*!< age of node */
516 Calibration *calibration; /*!< pointer to calibration data */
520 /* typedef for binary tree */
523 char name[100]; /*!< name of tree */
524 int memNodes; /*!< number of allocated nodes (do not exceed!) */
525 int nNodes; /*!< number of nodes in tree (including lower root in rooted trees) */
526 int nIntNodes; /*!< number of interior nodes in tree (excluding lower root in rooted trees) */
527 int isRooted; /*!< is tree rooted? */
528 int isClock; /*!< is tree clock? */
529 int isCalibrated; /*!< is tree calibrated? */
530 int nRelParts; /*!< number of relevant partitions */
531 int *relParts; /*!< pointer to relevant partitions */
532 int checkConstraints; /*!< does tree have constraints? */
533 int nConstraints; /*!< number of constraints */
534 int *constraints; /*!< pointer to constraints */
535 int nLocks; /*!< number of constrained (locked) nodes */
536 TreeNode **allDownPass; /*!< downpass array of all nodes */
537 TreeNode **intDownPass; /*!< downpass array of interior nodes (including upper but excluding lower root in rooted trees) */
538 TreeNode *root; /*!< pointer to root (lower root in rooted trees) */
539 TreeNode *nodes; /*!< array containing the nodes */
540 BitsLong *bitsets; /*!< pointer to bitsets describing splits */
541 BitsLong *flags; /*!< pointer to cond like flags */
542 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 */
546 /* typedef for node in polytomous tree */
549 char label[100]; /*!< name of node if terminal */
550 struct pNode *left, *sib, *anc; /*!< pointers to adjacent nodes */
551 int x, y, mark; /*!< scratch variables */
552 int partitionIndex; /*!< partition index in sumt (scratch) */
553 int index; /*!< index of node (if < numLocalTaxa = local taxon index) */
554 int memoryIndex; /*!< immutable index of memory position */
555 int isLocked; /*!< is the node locked? */
556 int lockID; /*!< id of lock */
557 int isDated; /*!< is node dated? */
558 MrBFlt length; /*!< age of node */
559 MrBFlt depth; /*!< depth (height) of node */
560 MrBFlt age; /*!< age of node */
561 MrBFlt support, f; /*!< scratch variables */
562 BitsLong *partition; /*!< pointer to partition (split) bitset */
563 Calibration *calibration; /*!< pointer to dating of node */
567 /* typedef for polytomous tree */
570 char name[100]; /*!< name of tree */
571 int memNodes; /*!< number of allocated nodes; do not exceed! */
572 int nNodes; /*!< number of nodes in tree */
573 int nIntNodes; /*!< number of interior nodes in tree */
574 PolyNode **allDownPass; /*!< downpass array over all nodes */
575 PolyNode **intDownPass; /*!< downpass array over interior nodes */
576 PolyNode *root; /*!< pointer to root (lower for rooted trees */
577 PolyNode *nodes; /*!< array holding the tree nodes */
578 BitsLong *bitsets; /*!< bits describing partitions (splits) */
579 int nBSets; /*!< number of effective branch length sets */
580 int nESets; /*!< number of breakpoint rate sets */
581 char **bSetName; /*!< names of effective branch length sets */
582 char **eSetName; /*!< names of breakpoint rate sets */
583 int **nEvents; /*!< number of branch events of bp rate set */
584 MrBFlt ***position; /*!< position of branch events */
585 MrBFlt ***rateMult; /*!< parameter of branch events */
586 MrBFlt **effectiveBrLen; /*!< effective branch lengths of ebl set */
587 int brlensDef; /*!< are brlens defined ? */
588 int isRooted; /*!< is tree rooted? */
589 int isClock; /*!< is tree clock? */
590 int isCalibrated; /*!< is tree calibrated? */
591 int isRelaxed; /*!< is tree relaxed? */
592 MrBFlt clockRate; /*!< clock rate */
593 int popSizeSet; /*!< does tree have a population size set? */
594 MrBFlt *popSize; /*!< the population size */
595 char *popSizeSetName; /*!< name of the population size set */
599 /* struct for holding model parameter info for the mcmc run */
602 int index; /* index to the parameter (0, 1, 2, ...) */
603 int paramType; /* the type of the parameter */
604 int paramId; /* unique ID for parameter x prior combination */
605 MrBFlt *values; /* main values of parameter */
606 MrBFlt *subValues; /* subvalues of parameter */
607 int *intValues; /* integer values (model index/growth fxn) */
608 int nValues; /* number of values */
609 int nSubValues; /* number of subvalues */
610 int nIntValues; /* number of intvalues */
611 MrBFlt min; /* minimum value of parameter */
612 MrBFlt max; /* maximum value of parameter */
613 int *relParts; /* pointer to relevant divisions */
614 int nRelParts; /* number of relevant divisions */
615 int upDate; /* update flag (for copying) */
616 struct param **subParams; /* pointers to subparams (for topology) */
617 int nSubParams; /* number of subparams */
618 Tree **tree; /* pointer to tree ptrs (for brlens & topology) */
619 int treeIndex; /* index to first tree in mcmcTree */
620 int hasBinaryStd; /* has binary standard chars */
621 int *sympiBsIndex; /* pointer to sympi bsIndex (std chars) */
622 int *sympinStates; /* pointer to sympi nStates (std chars) */
623 int *sympiCType; /* pointer to sympi cType (std chars) */
624 int nSympi; /* number of sympis */
625 int printParam; /* whether parameter should be printed */
626 int nPrintSubParams; /* number of subparams that should be printed */
627 char *paramHeader; /* a string holding header for param values */
628 char *name; /* string holding name of parameter */
629 char *paramTypeName; /* pointer to description of parameter type */
630 int checkConstraints; /* is tree parameter constrained? */
631 int fill; /* flags whether the parameter should be filled */
632 int nStdStateFreqs; /* number of std state frequencies */
633 MrBFlt *stdStateFreqs; /* pointer to std state frequencies */
634 int **nEvents; /* number of branch events for Cpp model */
635 /* nEvents[0..2*numCains][0..numNodes=2*numTaxa] */
636 MrBFlt ***position; /* event positions for Cpp relaxed clock model */
637 MrBFlt ***rateMult; /* rate multipliers for Cpp relaxed clock model */
638 int affectsLikelihood; /* does parameter directly influence likelihood? */
639 MrBFlt* priorParams; /* pointer to the prior parameters */
640 LnPriorProbFxn LnPriorProb; /* ln prior prob function */
641 LnPriorRatioFxn LnPriorRatio; /* ln prior prob ratio function */
644 #if defined(THREADS_ENABLED)
647 typedef struct s_launch_struct
655 /* parameter ID values */
656 /* identifies unique model parameter x prior combinations */
664 #define SYMPI_UNI_MS 8
666 #define SYMPI_EXP_MS 10
668 #define SYMPI_FIX_MS 12
669 #define SYMPI_EQUAL 13
672 #define PI_EMPIRICAL 16
678 #define PINVAR_UNI 22
679 #define PINVAR_FIX 23
680 #define CORREL_UNI 24
681 #define CORREL_FIX 25
682 #define SWITCH_UNI 26
683 #define SWITCH_EXP 27
684 #define SWITCH_FIX 28
685 #define RATEMULT_DIR 29
686 #define RATEMULT_FIX 30
687 #define TOPOLOGY_NCL_UNIFORM 31
688 #define TOPOLOGY_NCL_CONSTRAINED 32
689 #define TOPOLOGY_NCL_FIXED 33
690 #define TOPOLOGY_NCL_UNIFORM_HOMO 34
691 #define TOPOLOGY_NCL_UNIFORM_HETERO 35
692 #define TOPOLOGY_NCL_CONSTRAINED_HOMO 36
693 #define TOPOLOGY_NCL_CONSTRAINED_HETERO 37
694 #define TOPOLOGY_NCL_FIXED_HOMO 38
695 #define TOPOLOGY_NCL_FIXED_HETERO 39
696 #define TOPOLOGY_CL_UNIFORM 40
697 #define TOPOLOGY_CL_CONSTRAINED 41
698 #define TOPOLOGY_CL_FIXED 42
699 #define TOPOLOGY_CCL_UNIFORM 43
700 #define TOPOLOGY_CCL_CONSTRAINED 44
701 #define TOPOLOGY_CCL_FIXED 45
702 #define TOPOLOGY_PARSIMONY_UNIFORM 46
703 #define TOPOLOGY_PARSIMONY_CONSTRAINED 47
704 #define TOPOLOGY_PARSIMONY_FIXED 48
705 #define BRLENS_UNI 49
706 #define BRLENS_EXP 50
707 #define BRLENS_GamDir 51
708 #define BRLENS_iGmDir 52
709 #define BRLENS_twoExp 53
710 #define BRLENS_FIXED 54
711 #define BRLENS_CLOCK_UNI 55
712 #define BRLENS_CLOCK_COAL 56
713 #define BRLENS_CLOCK_BD 57
714 #define BRLENS_CLOCK_FIXED 58
715 #define BRLENS_CLOCK_SPCOAL 59
716 #define BRLENS_CLOCK_FOSSIL 60
717 #define BRLENS_PARSIMONY 61
718 #define SPECRATE_UNI 62
719 #define SPECRATE_EXP 63
720 #define SPECRATE_FIX 64
721 #define EXTRATE_BETA 65
722 #define EXTRATE_FIX 66
723 #define FOSLRATE_BETA 67
724 #define FOSLRATE_FIX 68
725 #define POPSIZE_UNI 69
726 #define POPSIZE_GAMMA 70
727 #define POPSIZE_FIX 71
728 #define POPSIZE_NORMAL 72
729 #define POPSIZE_LOGNORMAL 73
730 #define AAMODEL_FIX 74
731 #define AAMODEL_MIX 75
732 #define GROWTH_UNI 76
733 #define GROWTH_EXP 77
734 #define GROWTH_FIX 78
735 #define GROWTH_NORMAL 79
752 #define OMEGA_10UUB 96
753 #define OMEGA_10UUF 97
754 #define OMEGA_10UEB 98
755 #define OMEGA_10UEF 99
756 #define OMEGA_10UFB 100
757 #define OMEGA_10UFF 101
758 #define OMEGA_10EUB 102
759 #define OMEGA_10EUF 103
760 #define OMEGA_10EEB 104
761 #define OMEGA_10EEF 105
762 #define OMEGA_10EFB 106
763 #define OMEGA_10EFF 107
764 #define OMEGA_10FUB 108
765 #define OMEGA_10FUF 109
766 #define OMEGA_10FEB 110
767 #define OMEGA_10FEF 111
768 #define OMEGA_10FFB 112
769 #define OMEGA_10FFF 113
770 #define CPPRATE_FIX 114
771 #define CPPRATE_EXP 115
772 #define CPPMULTDEV_FIX 116
773 #define TK02VAR_FIX 117
774 #define TK02VAR_EXP 118
775 #define TK02VAR_UNI 119
776 #define TOPOLOGY_RCL_UNIFORM 120
777 #define TOPOLOGY_RCL_CONSTRAINED 121
778 #define TOPOLOGY_RCL_FIXED 122
779 #define TOPOLOGY_RCCL_UNIFORM 123
780 #define TOPOLOGY_RCCL_CONSTRAINED 124
781 #define TOPOLOGY_RCCL_FIXED 125
782 #define TOPOLOGY_SPECIESTREE 126
783 #define CPPEVENTS 127
784 #define TK02BRANCHRATES 128
785 #define TOPOLOGY_FIXED 129
786 #define IGRVAR_FIX 130
787 #define IGRVAR_EXP 131
788 #define IGRVAR_UNI 132
789 #define IGRBRANCHRATES 133
790 #define CLOCKRATE_FIX 134
791 #define CLOCKRATE_NORMAL 135
792 #define CLOCKRATE_LOGNORMAL 136
793 #define CLOCKRATE_GAMMA 137
794 #define CLOCKRATE_EXP 138
795 #define SPECIESTREE_UNIFORM 139
796 #define GENETREERATEMULT_DIR 140
797 #define GENETREERATEMULT_FIX 141
798 #define REVMAT_MIX 142
799 #define MIXEDVAR_FIX 143
800 #define MIXEDVAR_EXP 144
801 #define MIXEDVAR_UNI 145
802 #define MIXEDBRCHRATES 146
804 #if defined (BEAGLE_ENABLED)
805 #define MB_BEAGLE_SCALE_ALWAYS 0
806 #define MB_BEAGLE_SCALE_DYNAMIC 1
808 #define MB_PRINT_DYNAMIC_RESCALE_FAIL_STAT
812 /* typedef for a MoveFxn */
813 typedef int (MoveFxn)(Param *param, int chain, RandLong *seed, MrBFlt *lnPriorRatio, MrBFlt *lnProposalRatio, MrBFlt *mvp);
815 /* typedef for an ApplicFxn */
816 typedef int (ApplicFxn)(Param *param);
818 /* typedef for an AutotuneFxn */
819 typedef void (AutotuneFxn)(MrBFlt acceptanceRate, MrBFlt targetRate, int batch, MrBFlt *tuningParameter, MrBFlt minTuning, MrBFlt maxTuning);
821 /* struct holding info on each move type that the program handles */
824 MoveFxn *moveFxn; /* pointer to the move function */
825 ApplicFxn *isApplicable; /* pointer to function determining whether move is applicable to a parameter */
826 int nApplicable; /* number of relevant params */
827 int applicableTo[40]; /* pointer to ID of relevant params */
828 char *name; /* name of the move type */
829 char *shortName; /* abbreviated name of the move type */
830 char *paramName; /* name of subparameter if complex parameter */
831 int subParams; /* are we changing subparams (brlens of topol.) */
832 char *tuningName[5]; /* name of tuning params */
833 char *shortTuningName[5];/* short name of tuning params */
834 MrBFlt relProposalProb; /* default relative proposal probability */
835 int numTuningParams; /* number of tuning parameters */
836 MrBFlt tuningParam[5]; /* default tuning parameters for the proposal */
837 MrBFlt minimum[5]; /* minimum values for tuning params */
838 MrBFlt maximum[5]; /* maximum values for tuning params */
839 int parsimonyBased; /* this move is based on parsimony (YES/NO) */
840 int level; /* user level of this move */
841 AutotuneFxn *Autotune; /* pointer to the autotuning function */
842 MrBFlt targetRate; /* default target acceptance rate for autotuning */
845 /* max number of move types */
846 #define NUM_MOVE_TYPES 100
848 /* struct holding info on each move */
849 /* Note: This allows several different moves to affect the same parameter */
850 /* It also allows the same move to affect different parameters as before */
851 /* This is also a good place to keep the proposal probs */
854 char *name; /* name of the move */
855 MoveType *moveType; /* pointer to the move type */
856 MoveFxn *moveFxn; /* pointer to the move function */
857 Param *parm; /* ptr to parameter the move applies to */
858 MrBFlt *relProposalProb; /* the relative proposal probability */
859 MrBFlt *cumProposalProb; /* the cumulative proposal probability */
860 int *nAccepted; /* number of accepted moves */
861 int *nTried; /* number of tried moves */
862 int *nBatches; /* counter for autotuning rounds */
863 int *nTotAccepted; /* total number of accepted moves */
864 int *nTotTried; /* total number of tried moves */
865 MrBFlt *targetRate; /* target acceptance rate for autotuning */
866 MrBFlt *lastAcceptanceRate;/* acceptance rate in last complete batch */
867 MrBFlt **tuningParam; /* tuning parameters for the move */
870 typedef int (*LikeDownFxn)(TreeNode *, int, int);
871 typedef int (*LikeRootFxn)(TreeNode *, int, int);
872 typedef int (*LikeScalerFxn)(TreeNode *, int, int);
873 typedef int (*LikeFxn)(TreeNode *, int, int, MrBFlt *, int);
874 typedef int (*TiProbFxn)(TreeNode *, int, int);
875 typedef int (*LikeUpFxn)(TreeNode *, int, int);
876 typedef int (*PrintAncStFxn)(TreeNode *, int, int);
877 typedef int (*StateCodeFxn) (int);
878 typedef int (*PrintSiteRateFxn) (TreeNode *, int, int);
879 typedef int (*PosSelProbsFxn) (TreeNode *, int, int);
880 typedef int (*SiteOmegasFxn) (TreeNode *, int, int);
882 typedef struct cmdtyp
891 char *cmdDescription;
899 char *string; /* parameter name */
900 char *valueList; /* list of values that could be input */
901 ParmFxn fp; /* function pointer */
903 ParmInfo, *ParmInfoPtr;
907 int dataType; /* data type for partition */
908 int nStates; /* number of states for this type of data */
909 int codon[64]; /* gives protein ID for each codon */
910 int codonNucs[64][3]; /* gives triplet for each codon */
911 int codonAAs[64]; /* gives protein ID for implemented code */
913 char nucModel[100]; /* nucleotide model used */
914 char nst[100]; /* number of substitution types */
915 char parsModel[100]; /* use the (so-called) parsimony model */
916 char geneticCode[100]; /* genetic code used */
917 int coding; /* type of patterns encoded */
918 char codingString[100]; /* string describing type of patterns encoded */
919 char ploidy[100]; /* ploidy level */
920 char omegaVar[100]; /* type of omega variation model */
921 char ratesModel[100]; /* rates across sites model */
922 int numGammaCats; /* number of categories for gamma approximation */
923 char useGibbs[100]; /* flags whether Gibbs sampling of discrete gamma is used */
924 int gibbsFreq; /* frequency of Gibbs resampling of discrete gamma */
926 int numBetaCats; /* number of categories for beta approximation */
927 int numM10GammaCats; /* number of cats for gamma approx (M10 model) */
928 int numM10BetaCats; /* number of cats for beta approx (M10 model) */
929 char covarionModel[100];/* use covarion model? (yes/no) */
930 char augmentData[100]; /* should data be augmented */
932 char tRatioPr[100]; /* prior for ti/tv rate ratio */
935 char revMatPr[100]; /* prior for GTR model */
938 MrBFlt revMatSymDir; /* prior for mixed GTR subspace model */
939 char aaModelPr[100]; /* prior for amino acid model */
941 MrBFlt aaModelPrProbs[10];
942 char aaRevMatPr[100]; /* prior for aa GTR model */
943 MrBFlt aaRevMatFix[190];
944 MrBFlt aaRevMatDir[190];
945 char omegaPr[100]; /* prior for omega */
948 char ny98omega1pr[100]; /* prior for class 1 omega (Ny98 model) */
949 MrBFlt ny98omega1Fixed;
950 MrBFlt ny98omega1Beta[2];
951 char ny98omega3pr[100]; /* prior for class 3 omega (Ny98 model) */
952 MrBFlt ny98omega3Fixed;
953 MrBFlt ny98omega3Uni[2];
954 MrBFlt ny98omega3Exp;
955 char m3omegapr[100]; /* prior for all three omegas (M3 model) */
956 MrBFlt m3omegaFixed[3];
957 char m10betapr[100]; /* prior for omega variation (M10 model) */
958 char m10gammapr[100];
960 MrBFlt m10betaUni[2];
961 MrBFlt m10betaFix[2];
963 MrBFlt m10gammaUni[2];
964 MrBFlt m10gammaFix[2];
965 char codonCatFreqPr[100]; /* prior for selection cat frequencies */
966 MrBFlt codonCatFreqFix[3];
967 MrBFlt codonCatDir[3];
968 char stateFreqPr[100]; /* prior for character state frequencies */
969 MrBFlt stateFreqsFix[200];
970 MrBFlt stateFreqsDir[200];
971 char stateFreqsFixType[100];
973 char shapePr[100]; /* prior for gamma/lnorm shape parameter */
977 char pInvarPr[100]; /* prior for proportion of invariable sites */
980 char adGammaCorPr[100]; /* prior for correlation param of adGamma model */
983 char covSwitchPr[100]; /* prior for switching rates of covarion model */
984 MrBFlt covswitchFix[2];
985 MrBFlt covswitchUni[2];
987 char symPiPr[100]; /* prior for pi when unidentifiable states used */
989 MrBFlt symBetaUni[2];
991 char ratePr[100]; /* prior on rate for a partition */
993 char generatePr[100]; /* prior on rate for a gene (one or more partitions) */
994 MrBFlt generatePrDir;
995 char brownCorPr[100]; /* prior for correlation of Brownian model */
997 MrBFlt brownCorrUni[2];
998 char brownScalesPr[100]; /* prior for scales of Brownian model */
999 MrBFlt brownScalesFix;
1000 MrBFlt brownScalesUni[2];
1001 MrBFlt brownScalesGamma[2];
1002 MrBFlt brownScalesGammaMean;
1004 char topologyPr[100]; /* prior for tree topology */
1005 int topologyFix; /* user tree index for fixed topology */
1006 int *activeConstraints; /* which constraints are active? */
1007 int numActiveConstraints;
1009 char brlensPr[100]; /* prior on branch lengths */
1010 int brlensFix; /* user tree index for fixed brlens */
1011 MrBFlt brlensUni[2];
1013 MrBFlt brlens2Exp[2];
1014 MrBFlt brlensDir[4];
1015 MrBFlt brlensGamma[2];
1016 char speciesTreeBrlensPr[100]; /* prior on branch lengths of species tree */
1017 char unconstrainedPr[100]; /* prior on branch lengths if unconstrained */
1018 char clockPr[100]; /* prior on branch if clock enforced */
1019 char clockVarPr[100]; /* prior on clock rate variation (strict, cpp, tk02, igr, ...) */
1020 char nodeAgePr[100]; /* prior on node depths (unconstrained, constraints) */
1021 char speciationPr[100]; /* prior on speciation rate (net diversification) */
1022 MrBFlt speciationFix;
1023 MrBFlt speciationUni[2];
1024 MrBFlt speciationExp;
1025 char extinctionPr[100]; /* prior on relative extinction rate (turnover) */
1026 MrBFlt extinctionFix;
1027 MrBFlt extinctionBeta[2];
1028 char fossilizationPr[100]; /* prior on fossilization rate (sampling proportion) */
1029 MrBFlt fossilizationFix;
1030 MrBFlt fossilizationBeta[2];
1031 char sampleStrat[100]; /* taxon sampling strategy (for b-d process) */
1032 int sampleFSNum; /* number of fossil slice sampling events (s) */
1033 MrBFlt *sampleFSTime; /* fossil slice sampling times (t_i, i=1,..,s) */
1034 MrBFlt *sampleFSProb; /* fossil slice sampling probs (rho_i, i=1,..,s) */
1035 MrBFlt sampleProb; /* extant taxon sampling fraction (rho) */
1036 Calibration treeAgePr; /* prior on tree age for uniform clock trees */
1037 char clockRatePr[100]; /* prior on base substitution rate of tree for clock trees */
1038 MrBFlt clockRateNormal[2];
1039 MrBFlt clockRateLognormal[2];
1040 MrBFlt clockRateGamma[2];
1041 MrBFlt clockRateExp;
1042 MrBFlt clockRateFix;
1043 char popSizePr[100]; /* prior on population size */
1045 MrBFlt popSizeUni[2];
1046 MrBFlt popSizeLognormal[2];
1047 MrBFlt popSizeGamma[2];
1048 MrBFlt popSizeNormal[2];
1049 char popVarPr[100]; /* prior on pop. size variation across tree */
1050 char growthPr[100]; /* prior on coalescence growth rate */
1052 MrBFlt growthUni[2];
1054 MrBFlt growthNorm[2];
1055 char cppRatePr[100]; /* prior on CPP rate */
1058 char cppMultDevPr[100]; /* prior on CPP rate multiplier Lognormal variance */
1059 MrBFlt cppMultDevFix;
1060 char tk02varPr[100]; /* prior on TK02 lognormal rate variance */
1062 MrBFlt tk02varUni[2];
1064 char igrvarPr[100]; /* prior on IGR gamma distribution variance */
1066 MrBFlt igrvarUni[2];
1068 char mixedvarPr[100]; /* prior on mixed relaxed clock rate variance */
1070 MrBFlt mixedvarUni[2];
1073 char tratioFormat[30]; /* format used to report tratio */
1074 char revmatFormat[30]; /* format used to report revmat */
1075 char ratemultFormat[30]; /* format used to report ratemult */
1076 char treeFormat[30]; /* format used to report trees/topologies */
1077 char inferAncStates[5]; /* should ancestral states be inferred (Yes/No)? */
1078 char inferSiteOmegas[5]; /* should site omega vals be inferred (Yes/No)? */
1079 char inferSiteRates[5]; /* should site rates be inferred (Yes/No)? */
1080 char inferPosSel[5]; /* should site selection be inferred (Yes/No)? */
1081 } Model, ModelParams;
1083 typedef struct chain
1085 int numGen; /* number of MCMC cycles */
1086 int sampleFreq; /* frequency to sample chain */
1087 int printFreq; /* frequency to print chain */
1088 int swapFreq; /* frequency to attempt swap of states */
1089 int numRuns; /* number of runs */
1090 int numChains; /* number of chains */
1091 int isSS; /* do we do Steppingstone Sampling */
1092 int startFromPriorSS; /* If Yes SS is moving from Prior to Posterior */
1093 int numStepsSS; /* Number of steps in SS */
1094 int burninSS; /* Fixed burnin for SS */
1095 MrBFlt alphaSS; /* Beta values are distributed according to quantiles of Beta(alphaSS,1.0) distribution */
1096 int backupCheckSS; /* Frequency of checkpoints backup */
1097 MrBFlt chainTemp; /* chain temperature */
1098 int userDefinedTemps; /* should we use the users temperatures? */
1099 MrBFlt userTemps[MAX_CHAINS]; /* user-defined chain temperatures */
1100 char chainFileName[100]; /* chain file name for output */
1101 int chainBurnIn; /* chain burn in length */
1102 int numStartPerts; /* number of perturbations to starting tree */
1103 char startTree[100]; /* starting tree for chain (current/random) */
1104 char startParams[100]; /* starting values for chain (current/reset) */
1105 int saveBrlens; /* should branch lengths be saved */
1106 MrBFlt weightScheme[3]; /* percent chars to increase/decrease in weight */
1107 int calcPbf; /* should we calculate the pseudo Bayes factor */
1108 int pbfInitBurnin; /* initial burnin when calculating pseudo BF */
1109 int pbfSampleFreq; /* sample frequency for pseudo BF */
1110 int pbfSampleTime; /* how many cycles to calcualate site prob. */
1111 int pbfSampleBurnin; /* burnin period for each site for pseudo BF */
1112 int swapAdjacentOnly; /* whether we only swap adjacent temperatures */
1113 int redirect; /* should output go to stdout */
1114 int allChains; /* should stats be output for all chains? */
1115 int numSwaps; /* number of swaps to try each time */
1116 int mcmcDiagn; /* should mcmc diagnostics be output? */
1117 int diagnFreq; /* mcmc diagnostics frequency */
1118 int diagnStat; /* statistic to use for mcmc diagnostics */
1119 int relativeBurnin; /* should a relative burnin be used ? */
1120 MrBFlt burninFraction; /* the sample fraction to discard as burnin */
1121 int allComps; /* top conv diagnosis for all pairs? */
1122 MrBFlt minPartFreq; /* minimum partition frequency for conv diagn */
1123 MrBFlt stopVal; /* top conv diagn value to reach before stopping */
1124 int stopRule; /* use stop rule? */
1125 STATS *stat; /* ptr to structs with mcmc diagnostics info */
1126 Tree *dtree; /* pointing to tree used for conv diagnostics */
1127 TreeList *treeList; /* vector of tree lists for saving trees */
1128 int saveTrees; /* save tree samples for later removal? */
1129 int stopTreeGen; /* generation after which no trees need be saved */
1130 fpos_t *tFilePos; /* position for reading trees for removal */
1131 int printMax; /* maximum number of chains to print */
1132 int printAll; /* whether to print all or only cold chains */
1133 int checkPoint; /* should we use check-pointing? */
1134 int checkFreq; /* check-pointing frequency */
1135 int runWithData; /* should we run with or without data? */
1136 int orderTaxa; /* order taxa before printing tree to file? */
1137 int append; /* order taxa before printing tree to file? */
1138 int autotune; /* autotune tuning parameters of proposals ? */
1139 int tuneFreq; /* autotuning frequency */
1142 typedef struct modelinfo
1144 /* General model information */
1145 int dataType; /* data type for partition */
1146 int nucModelId; /* structure of nucleotide model */
1147 int nst; /* # substitution types */
1148 int aaModelId; /* amino acid model type */
1149 int parsModelId; /* is parsimony model used YES/NO */
1151 /* Specific model information */
1152 int numGammaCats; /* number of gamma cats (1 if inapplic.) */
1153 int numBetaCats; /* number of beta cats (1 if inapplic.) */
1154 int numOmegaCats; /* number of omega cats (1 if inapplic.) */
1155 int numTiCats; /* number of cats needing different tis */
1156 int numModelStates; /* number of states including hidden ones */
1157 int numStates; /* number of observable discrete states */
1159 /* Model parameter pointers */
1160 Param *tRatio; /* ptr to tRatio used in model */
1161 Param *revMat; /* ptr to revMat used in model */
1162 Param *omega; /* ptr to omega used in model */
1163 Param *stateFreq; /* ptr to statFreq used in model */
1164 Param *shape; /* ptr to shape used in model */
1165 Param *pInvar; /* ptr to pInvar used in model */
1166 Param *correlation; /* ptr to correlation used in model */
1167 Param *switchRates; /* ptr to switchRates (off<->on) */
1168 Param *rateMult; /* ptr to parition rateMult used in model */
1169 Param *geneTreeRateMult; /* ptr to gene tree rateMult used in model */
1170 Param *speciationRates; /* ptr to speciationRates used in model */
1171 Param *extinctionRates; /* ptr to extinctionRates used in model */
1172 Param *fossilizationRates; /* ptr to fossilizationRates */
1173 Param *popSize; /* ptr to population size used in model */
1174 Param *growthRate; /* ptr to growth rate used in model */
1175 Param *topology; /* ptr to topology used in model */
1176 Param *brlens; /* ptr to brlens (and tree) used in model */
1177 Param *speciesTree; /* ptr to species tree used in model */
1178 Param *aaModel; /* ptr to amino acid matrix used */
1179 Param *cppMultDev; /* ptr to cpp ratemult lognormal variance */
1180 Param *cppRate; /* ptr to CPP rate used in model */
1181 Param *cppEvents; /* ptr to CPP events */
1182 Param *tk02var; /* ptr to variance for TK02 relaxed clock */
1183 Param *tk02BranchRates; /* ptr to branch rates for TK02 relaxed clock */
1184 Param *igrvar; /* ptr to gamma var for IGR relaxed clock */
1185 Param *igrBranchRates; /* ptr to branch rates for IGR relaxed clock*/
1186 Param *mixedvar; /* ptr to var for mixed relaxed clock */
1187 Param *mixedBrchRates; /* ptr to branch rates for mixed relaxed clock */
1188 Param *clockRate; /* ptr to clock rate parameter */
1190 /* Information about characters and transformations */
1191 int numChars; /* number of compressed characters */
1192 int numUncompressedChars; /* number of uncompressed characters */
1193 int numDummyChars; /* number of dummy characters */
1194 int compMatrixStart; /* start column in compressed matrix */
1195 int compMatrixStop; /* stop column in compressed matrix */
1196 int compCharStart; /* start char among compressed chars */
1197 int compCharStop; /* stop char among compressed chars */
1198 int parsMatrixStart; /* start column in parsimony matrix */
1199 int parsMatrixStop; /* stop collumn in parsimony matrix */
1200 int nParsIntsPerSite; /* # parsimony ints per character */
1201 int nCharsPerSite; /* number chars per site (eg 3 for codon) */
1202 int rateProbStart; /* start of rate probs (for adgamma) */
1204 /* Variables for eigen decompositions */
1205 int cijkLength; /* stores length of cijk vector */
1206 int nCijkParts; /* stores number of cijk partitions (1 except for omega/covarion models) */
1207 int upDateCijk; /* whether cijk vector needs to be updated */
1209 /* Variables for standard model */
1210 int *tiIndex; /* index to trans probs for each compressed char*/
1211 int *bsIndex; /* index to stat freqs for each compressed char */
1212 int *nStates; /* # states of each compressed char */
1213 int *cType; /* whether char is ord, unord or irrev */
1214 int *weight; /* prior weight of each compressed char */
1215 int isTiNeeded[20]; /* marks whether a trans prob matrix is needed */
1217 /* Gibbs sampling of gamma site rate parameters */
1218 CLFlt ***catLike; /* likelihood for Gibbs sampling of gamma */
1219 CLFlt ***catLnScaler; /* scaler for Gibbs sampling of gamma */
1220 int gibbsGamma; /* flags whether Gibbs sampling of discrete gamma is used */
1221 int gibbsFreq; /* frequency of Gibbs resampling of discrete gamma */
1223 /* Variables for parsimony sets and parsimony calculations */
1224 MrBFlt parsTreeLength[MAX_CHAINS*2]; /* parsimony tree lengths for chains */
1225 BitsLong **parsSets; /* parsimony sets */
1226 int numParsSets; /* number of parsimony sets */
1227 CLFlt *parsNodeLens; /* parsimony node lengths */
1228 int numParsNodeLens; /* number of parsimony node lengths */
1230 /* Miscellaneous parameters */
1231 int mark; /* scratch parameter */
1232 int parsimonyBasedMove; /* is parsimony-based move used (YES/NO) */
1234 /* Variables for conditional likelihoods */
1235 int upDateCl; /* flags whether update of cond likes needed */
1236 int upDateAll; /* flags whether update of entire tree is needed*/
1237 int *isPartAmbig; /* is tip partially ambiguous? */
1238 int **termState; /* index arrays for terminal branch shortcuts */
1239 CLFlt *invCondLikes; /* space for the invariable cond likes */
1240 CLFlt **condLikes; /* space for the cond likes */
1241 CLFlt **tiProbs; /* space for the ti probs */
1242 CLFlt **scalers; /* space for the node and site scalers */
1243 CLFlt **clP; /* handy pointers to cond likes for ti cats */
1244 #if defined (SSE_ENABLED)
1245 __m128 **clP_SSE; /* handy pointers to cond likes, SSE version */
1246 int numSSEChars; /* number of compact SSE character groups */
1247 CLFlt *lnL_SSE; /* temp storage for log site likes */
1248 CLFlt *lnLI_SSE; /* temp storage for log site invariable likes */
1250 MrBFlt **cijks; /* space for cijks */
1251 int **condLikeIndex; /* index to cond like space for nodes & chains */
1252 int *condLikeScratchIndex; /* index to scratch space for node cond likes */
1253 int **nodeScalerIndex; /* index to scaler space for nodes & chains */
1254 int *nodeScalerScratchIndex; /* index to scratch space for node scalers */
1255 int **scalersSet; /* flags whether scalers are set */
1256 int *scalersSetScratch; /* scratch flag for whether scalers are set */
1257 int *siteScalerIndex; /* index to site scaler space for chains */
1258 int siteScalerScratchIndex; /* index to scratch space for site scalers */
1259 int **tiProbsIndex; /* index to ti probs for branches & chains */
1260 int *tiProbsScratchIndex; /* index to scratch space for branch ti probs */
1261 int *cijkIndex; /* index to cijks for chains */
1262 int cijkScratchIndex; /* index to scratch space for cijks */
1263 int numCondLikes; /* number of cond like arrays */
1264 int numScalers; /* number of scaler arrays */
1265 int numTiProbs; /* number of ti prob arrays */
1266 int condLikeLength; /* length of cond like array (incl. ti cats) */
1267 int tiProbLength; /* length of ti prob array */
1268 MrBFlt lnLike[MAX_CHAINS]; /* log like for chain */
1269 CLFlt *ancStateCondLikes; /* ancestral state cond like array */
1271 /* Likelihood function pointers */
1272 LikeDownFxn CondLikeDown; /* function for calculating partials */
1273 LikeRootFxn CondLikeRoot; /* function for calculating partials at root */
1274 LikeScalerFxn CondLikeScaler; /* function for scaling partials */
1275 LikeFxn Likelihood; /* function for getting cond likes for tree */
1276 TiProbFxn TiProbs; /* function for calculating transition probs */
1277 LikeUpFxn CondLikeUp; /* final-pass calculation of cond likes */
1278 PrintAncStFxn PrintAncStates; /* function for sampling ancestral states */
1279 StateCodeFxn StateCode; /* function for getting states from codes */
1280 PrintSiteRateFxn PrintSiteRates; /* function for samling site rates */
1281 PosSelProbsFxn PosSelProbs; /* function for sampling pos. selection probs */
1282 SiteOmegasFxn SiteOmegas; /* function for sampling site omega values */
1284 /* Report variables */
1285 int printAncStates; /* should ancestral states be printed (YES/NO) */
1286 int printSiteRates; /* should site rates be printed (YES/NO) */
1287 int printPosSel; /* should selection be printed (YES/NO) */
1288 int printSiteOmegas; /* should site omegas be printed (YES/NO) */
1290 /* likelihood calculator flags */
1291 int useBeagle; /* use Beagle for this partition? */
1292 int useSSE; /* use SSE for this partition? */
1294 #if defined (BEAGLE_ENABLED)
1295 /* Beagle variables */
1296 int useBeagleResource; /* try to use this BEAGLE resource number */
1297 MrBFlt* branchLengths; /* array of branch lengths for Beagle */
1298 MrBFlt* inRates; /* array of category rates for Beagle */
1299 int* tiProbIndices; /* array of trans prob indices for Beagle */
1300 MrBFlt* logLikelihoods; /* array of log likelihoods from Beagle */
1301 int beagleInstance; /* beagle instance for division */
1302 MrBFlt* inWeights; /* array of weights for Beagle root likelihood */
1303 int* bufferIndices; /* array of partial indices for root likelihood */
1304 int* eigenIndices; /* array of eigen indices for root likelihood */
1305 int* childBufferIndices; /* array of child partial indices (unrooted) */
1306 int* childTiProbIndices; /* array of child ti prob indices (unrooted) */
1307 int* cumulativeScaleIndices; /* array of cumulative scale indices */
1308 int rescaleBeagleAll; /* set to rescale all nodes */
1309 int* rescaleFreq; /* rescale frequency for each chain's tree */
1310 int rescaleFreqOld; /* holds rescale frequency of current state */
1311 int recalculateScalers; /* shoud we recalculate scalers for current state YES/NO */
1312 int* succesCount; /* count number of succesful computation since last reset of scalers */
1313 int** isScalerNode; /* for each node and chain set to YES if scaled node */
1314 int* isScalerNodeScratch; /* scratch space to hold isScalerNode of proposed state*/
1315 long* beagleComputeCount; /* count of number of calls to likelihood */
1322 int *absentTaxa; /* information on absent taxa */
1323 int brlensDef; /* branch lengths defined? */
1324 char sumtFileName[100]; /* name of input file */
1325 char sumtOutfile[120]; /* name of output file */
1326 char curFileName[120]; /* name of file being processed */
1327 int burnin; /* actual burnin when parsing tree files */
1328 char sumtConType[100]; /* consensus tree type */
1329 int calcTreeprobs; /* should the individual tree probs be calculated*/
1330 int showSumtTrees; /* should the individual tree probs be shown */
1331 int printBrlensToFile; /* should branch lengths be printed to file */
1332 MrBFlt brlensFreqDisplay; /* threshold for printing branch lengths to file */
1333 int numRuns; /* number of independent analyses to summarize */
1334 int numTrees; /* number of tree params to summarize */
1335 int orderTaxa; /* order taxa in trees? */
1336 MrBFlt minPartFreq; /* minimum part. freq. for overall diagnostics */
1337 int table; /* show table of partition frequencies? */
1338 int summary; /* show summary diagnostics ? */
1339 int showConsensus; /* show consensus trees ? */
1340 int consensusFormat; /* format of consensus tree */
1341 PolyTree *tree; /* for storing tree read from file */
1342 int *order; /* for storing topology read from file */
1343 int orderLen; /* length of order array */
1344 int numTreesInLastBlock; /* number of trees in last block */
1345 int numTreesEncountered; /* number of trees encounted in total */
1346 int numTreesSampled; /* number of sampled trees in total */
1347 int isRooted; /* is sumt tree rooted ? */
1348 int isRelaxed; /* is sumt tree a relaxed clock tree ? */
1349 int isClock; /* is sumt tree a clock tree ? */
1350 int isCalibrated; /* is sumt tree calibrated ? */
1351 int nESets; /* number of event sets */
1352 int nBSets; /* number of branch rate sets */
1353 char **bSetName; /* name of effective branch length sets */
1354 char **eSetName; /* name of event sets */
1355 int popSizeSet; /* do sumt trees have population size set? */
1356 char *popSizeSetName; /* name of population size set */
1357 int BitsLongsNeeded; /* number of safe longs needed for taxon bits */
1358 int runId; /* id of run being processed */
1359 int numTaxa; /* number of sumt taxa */
1360 char **taxaNames; /* names of sumt taxa */
1361 int *numFileTrees; /* number of trees per file */
1362 int *numFileTreesSampled; /* number of trees sampled per file */
1363 int HPD; /* use highest posterior density? */
1366 /* formats for consensus trees */
1370 typedef struct comptree
1372 char comptFileName1[120]; /* name of first input file */
1373 char comptFileName2[120]; /* name of second input file */
1374 char comptOutfile[120]; /* name of output file */
1375 int burnin; /* actual burnin used when parsing tree files */
1376 MrBFlt minPartFreq; /* use partitions with frequency >= minPartFreq */
1381 char sumpFileName[100]; /* name of input file */
1382 char sumpOutfile[120]; /* name of output file */
1383 //int plot; /* output plot (y/n)? */
1384 int table; /* output table (y/n)? */
1385 int margLike; /* output marginal likelihood (y/n)? */
1386 int numRuns; /* number of independent analyses to summarize */
1387 int allRuns; /* should data for all runs be printed (yes/no)? */
1388 int HPD; /* use highest posterior density? */
1389 MrBFlt minProb; /* cut-off for model probabilities to show */
1392 typedef struct sumss
1394 //int plot; /* output plot (y/n)? */
1395 int numRuns; /* number of independent analyses to summarize */
1396 int allRuns; /* should data for all runs be printed (yes/no)? */
1397 int stepToPlot; /* Which step to plot in the step plot */
1398 int askForMorePlots; /* Should user be asked to plot for different discardfraction (y/n)? */
1399 int smoothing; /* An integer indicating number of neighbors to average over when dooing smoothing of curvs on plots */
1400 MrBFlt discardFraction; /* Proportion of samples discarded when ploting step plot.*/
1405 char plotFileName[120]; /* name of input file */
1406 char parameter[100]; /* parameter(s) to be plotted */
1407 char match[100]; /* whether the match needs to be perfect */
1412 int numTrees; /* number of trees to reassemble */
1413 int numRuns; /* number of runs to reassemble */
1416 typedef struct doublet
1418 BitsLong first, second;
1421 typedef struct matrix
1430 typedef struct charinfo
1443 int isExcluded; /* is the character excluded */
1444 int numStates; /* number of observed states for the character */
1445 int charType; /* type of character */
1446 int isMissAmbig; /* is the character missing or ambiguous */
1447 int ctype; /* ordering of character */
1448 int charId; /* char ID index for doublet and codon models */
1449 int pairsId; /* char ID for doublets */
1450 int bigBreakAfter; /* is there a large break after this character */
1456 int isDeleted; /* is the taxon deleted */
1457 int charCount; /* count holder */
1486 /* global variables */
1487 extern int abortMove; /* flag determining whether to abort move */
1488 extern int *activeParams[NUM_LINKED]; /* a table holding the parameter status */
1489 extern int *activeParts; /* partitions changes should apply to */
1490 extern int autoClose; /* autoclose */
1491 extern int autoOverwrite; /* Overwrite or append outputfiles when nowarnings=yes */
1492 extern int chainHasAdgamma; /* indicates if chain has adgamma HMMs */
1493 extern Chain chainParams; /* holds parameters for Markov chain */
1494 extern CharInformation *charInfo; /* holds critical information about characters */
1495 extern char **charSetNames; /* holds names of character sets */
1496 extern int *compCharPos; /* char position in compressed matrix */
1497 extern int *compColPos; /* column position in compressed matrix */
1498 extern BitsLong *compMatrix; /* compressed character matrix */
1499 extern int compMatrixRowSize; /* row size of compressed matrix */
1500 extern Comptree comptreeParams; /* holds parameters for comparetree command */
1501 extern char **constraintNames; /* holds names of constraints */
1502 extern BitsLong **definedConstraint; /* holds information about defined constraints */
1503 extern BitsLong **definedConstraintTwo; /* bitfields representing second taxa sets of defined constraints (for PARTIAL constraints) */
1504 extern BitsLong **definedConstraintPruned; /* bitfields representing taxa sets of defined constraints after delited taxa are removed */
1505 extern BitsLong **definedConstraintTwoPruned; /* bitfields representing second taxa sets of defined constraints after delited taxa are removed (for PARTIAL constraints) */
1506 extern int dataType; /* type of data */
1507 extern Calibration defaultCalibration; /* default model settings */
1508 extern ModelParams defaultModel; /* default model settings */
1509 extern int defChars; /* flag for whether number of characters is known*/
1510 extern int defMatrix; /* flag for whether matrix is successfull read */
1511 extern int defPairs; /* flag for whether constraints on tree are read */
1512 extern int defPartition; /* flag for whether character partition is read */
1513 extern int defTaxa; /* are taxon labels defined ? */
1514 extern Doublet doublet[16]; /* holds information on states for doublets */
1515 extern int echoMB; /* flag used by Manual to prevent echoing */
1516 extern BitsLong expecting; /* variable denoting expected token type */
1517 extern int fileNameChanged; /* has file name been changed? */
1518 extern int foundNewLine; /* whether a new line has been found */
1519 extern char gapId; /* gap character Id */
1520 extern RandLong globalSeed; /* seed that is initialized at start up */
1521 extern char **headerNames; /* string to hold headers in sump and plot */
1522 extern int inComment; /* flag for whether input stream is commented */
1523 extern int inferAncStates; /* should ancestral states be inferred (y/n) */
1524 extern int inferSiteOmegas; /* should site omega values be inferred (y/n) */
1525 extern int inferSiteRates; /* should site rates be inferred (y/n) */
1526 extern int inferPosSel; /* should positive selection be inferred (y/n) */
1527 extern char inputFileName[100]; /* input (NEXUS) file name */
1528 extern int inTreesBlock; /* are we in the sumt block */
1529 extern int inValidCommand; /* a useful flag set whenever you enter a cmd */
1530 extern int isInAmbig, isInPoly; /* flags whether we are within () or {} */
1531 extern int isMixed; /* flags whether dataset is mixed */
1532 extern int inMrbayesBlock; /* flag for whether we are in a mrbayes block */
1533 extern int *intValues; /* integer values of parameters */
1534 extern int isTaxsetDef; /* is a taxon set defined */
1535 extern int isTranslateDef; /* is a translation block defined */
1536 extern int isTranslateDiff; /* is translate different from current taxaset? */
1537 extern int *linkTable[NUM_LINKED]; /* how parameters are linked across parts */
1538 extern int localOutGroup; /* outgroup for non-excluded taxa */
1539 extern char **localTaxonNames; /* points to names of non-excluded taxa */
1540 extern FILE *logFileFp; /* file pointer to log file */
1541 extern char logFileName[100]; /* name of the log file */
1542 extern int logToFile; /* should screen output be logged to a file */
1543 extern char manFileName[100]; /* name of man file */
1544 extern char matchId; /* mach character Id */
1545 extern int *matrix; /* matrix containing original data */
1546 extern int matrixHasPoly; /* flag for whether matrix has polymorphisms */
1547 extern int memAllocs[NUM_ALLOCS]; /* allocated memory flags */
1548 extern int mode; /* mode of program (interactive/noninteractive) */
1549 extern char **modelIndicatorParams; /* model indicator params */
1550 extern char ***modelElementNames; /* names for component models */
1551 extern MCMCMove **moves; /* vector of applicable moves */
1552 extern MoveType moveTypes[NUM_MOVE_TYPES]; /* holds information on the move types */
1553 extern char missingId; /* missing character Id */
1554 extern Tree **mcmcTree; /* pointers to mcmc trees */
1555 extern Model *modelParams; /* holds model params for partitions */
1556 extern ModelInfo *modelSettings; /* stores important info on model params */
1557 extern int nBitsInALong; /* number of bits in a BitsLong */
1558 extern Calibration *nodeCalibration; /* holds information about node calibrations */
1559 extern int noWarn; /* no warnings on overwriting files */
1560 extern int nPThreads; /* number of pthreads to use */
1561 extern int numActiveLocks; /* number of active, locked nodes */
1562 extern int numApplicableMoves; /* number of moves applicable to parameters */
1563 extern int numChar; /* number of characters in character matrix */
1564 extern int numCharSets; /* holds number of character sets */
1565 extern int numComments; /* number of nested comments */
1566 extern int numCompressedChars; /* number of compressed characters */
1567 extern int numCurrentDivisions; /* number of partitions of data */
1568 extern int numDefinedConstraints; /* number of constraints defined */
1569 extern enum ConstraintType *definedConstraintsType; /* Store type of constraint */
1570 extern int numDefinedPartitions; /* number of partitions defined */
1571 extern int numDefinedSpeciespartitions; /* number of species partitions defined */
1572 extern int numGlobalChains; /* number of global chains */
1573 extern int numLocalTaxa; /* number of non-excluded taxa */
1574 extern int numLocalChar; /* number of non-excluded characters */
1575 extern int numMoveTypes; /* the number of move types */
1576 extern int numOpenExeFiles; /* number of execute files open */
1577 extern int numParams; /* number of parameters in model */
1578 extern int numDivisions; /* number of current divisions */
1579 extern int numPrintParams; /* number of substitution model parameters to print */
1580 extern int numPrintTreeParams; /* number of tree model parameters to print */
1581 extern CLFlt *numSitesOfPat; /* no. sites of each pattern */
1582 extern int numSpecies; /* number of species in current speciespartition */
1583 extern int numTaxa; /* number of taxa in character matrix */
1584 extern int numTaxaSets; /* holds number of taxa sets */
1585 extern int numTopologies; /* number of topologies for one chain and state */
1586 extern int numTranslates; /* number of taxa in active translate block */
1587 extern int numTrees; /* number of trees for one chain and state */
1588 extern int numUserTrees; /* number of defined user trees */
1589 extern int *numVars; /* number of variables in setting arrays */
1590 extern int *origChar; /* index from compressed char to original char */
1591 extern int outGroupNum; /* number of outgroup taxon */
1592 extern ParmInfo paramTable[]; /* information on parameters */
1593 extern MrBFlt *paramValues; /* values of parameters */
1594 extern int **partitionId; /* holds information about defined partitions */
1595 extern char **partitionNames; /* hold names of partitions (first is "default") */
1596 extern MrBFlt *parameterValues; /* vector holding sump or plot parameters */
1597 extern Param *params; /* vector of parameters in model */
1598 extern int partitionNum; /* index of current partition */
1599 extern Plot plotParams; /* holds parameters for plot command */
1600 extern int precision; /* precision of samples and summary stats */
1601 extern int *printAncStates; /* divisions to print anc states for */
1602 extern int quitOnError; /* quit on error? */
1603 extern int readComment; /* should we read comment (looking for &)? */
1604 extern int readWord; /* should we read a word next? */
1605 extern ReassembleInfo reassembleParams; /* holds parameters for reassemble command */
1606 extern int replaceLogFile; /* should logfile be replace/appended to */
1607 extern RandLong runIDSeed; /* seed used only for generating run ID [stamp] */
1608 extern BitsLong bitsLongWithAllBitsSet; /* a BitsLong with all bits set, for bit ops */
1609 extern int setUpAnalysisSuccess; /* Set to YES if analysis is set without error */
1610 extern int scientific; /* use scientific format for samples ? */
1611 extern ShowmovesParams showmovesParams; /* holds parameters for Showmoves command */
1612 extern char spacer[10]; /* holds blanks for printing indentations */
1613 extern NameSet *speciesNameSets; /* hold species name sets, one for each speciespartition */
1614 extern int **speciespartitionId; /* holds info about defined speciespartitions */
1615 extern char **speciespartitionNames; /* hold names of speciespartitions (first is "default") */
1616 extern int speciespartitionNum; /* index of current species partition */
1617 extern char stamp[11]; /* holds a unique identifier for each analysis */
1618 extern RandLong swapSeed; /* seed used only for determining which to swap */
1619 extern int state[MAX_CHAINS]; /* state of chain */
1620 extern MrBFlt *stdStateFreqs; /* std char state frequencies */
1621 extern int *stdType; /* compressed std char type: ord, unord, irrev */
1622 extern Sump sumpParams; /* holds parameters for sump command */
1623 extern char sumpToken[]; /* string holding a .p file token */
1624 extern char *sumpTokenP; /* pointer to a .p file token */
1625 extern Sumt sumtParams; /* holds parameters for sumt command */
1626 extern Sumss sumssParams; /* holds parameters for sumss command */
1627 extern int stdStateFreqsRowSize; /* row size for stdStateFreqs */
1628 extern int *sympiIndex; /* sympi state freq index for multistate chars */
1629 extern TaxaInformation *taxaInfo; /* holds critical information about taxa */
1630 extern char **taxaNames; /* holds name of taxa */
1631 extern char **taxaSetNames; /* holds names of taxa sets */
1632 extern BitsLong **taxaSet; /* holds information about defined taxasets */
1633 extern int *tempActiveConstraints; /* info on the active constraints in prset */
1634 extern int *tempLinkUnlink[NUM_LINKED]; /* for changing parameter linkage */
1635 extern int *tempLinkUnlinkVec; /* for changing parameter linkage */
1636 extern MrBFlt *tempNum; /* vector of numbers used for setting arrays */
1637 extern int *tempSet; /* temporarily holds defined character set */
1638 extern int theAmbigChar; /* int containing ambiguous character */
1639 extern int *tiIndex; /* compressed std char ti index */
1640 extern Calibration *tipCalibration; /* holds tip calibrations */
1641 extern char **transFrom; /* translation block information */
1642 extern char **transTo; /* translation block information */
1643 extern int userBrlensDef; /* are the branch lengths on user tree defined */
1644 extern int userLevel; /* the level of the user */
1645 extern PolyTree *userTree[]; /* array of user trees */
1646 extern char workingDir[100]; /* working directory */
1647 #if defined (BEAGLE_ENABLED)
1648 extern int tryToUseBEAGLE; /* try to use the BEAGLE library */
1649 extern long beagleFlags; /* BEAGLE requirement flags */
1650 extern int beagleResourceNumber; /* BEAGLE resource number */
1651 extern int* beagleResource; /* BEAGLE resource list */
1652 extern int beagleResourceCount; /* BEAGLE resource list length */
1653 extern int beagleInstanceCount; /* total number of BEAGLE instances */
1654 extern int beagleScalingScheme; /* BEAGLE dynamic scaling */
1655 extern int beagleScalingFrequency; /* BEAGLE rescaling frequency */
1656 extern int recalcScalers; /* shoud we recalculate scalers for one of divisions for current state YES/NO */
1658 #if defined (THREADS_ENABLED)
1659 extern int tryToUseThreads; /* try to use pthreads with BEAGLE library */
1662 /* Aamodel parameters */
1663 extern MrBFlt aaJones[20][20]; /* rates for Jones model */
1664 extern MrBFlt aaDayhoff[20][20]; /* rates for Dayhoff model */
1665 extern MrBFlt aaMtrev24[20][20]; /* rates for mtrev24 model */
1666 extern MrBFlt aaMtmam[20][20]; /* rates for mtmam model */
1667 extern MrBFlt aartREV[20][20]; /* rates for rtREV model */
1668 extern MrBFlt aaWAG[20][20]; /* rates for WAG model */
1669 extern MrBFlt aacpREV[20][20]; /* rates for aacpREV model */
1670 extern MrBFlt aaVt[20][20]; /* rates for VT model */
1671 extern MrBFlt aaBlosum[20][20]; /* rates for Blosum62 model */
1672 extern MrBFlt aaLG[20][20]; /* rates for LG model */
1673 extern MrBFlt jonesPi[20]; /* stationary frequencies for Jones model */
1674 extern MrBFlt dayhoffPi[20]; /* stationary frequencies for Dayhoff model */
1675 extern MrBFlt mtrev24Pi[20]; /* stationary frequencies for mtrev24 model */
1676 extern MrBFlt mtmamPi[20]; /* stationary frequencies for mtmam model */
1677 extern MrBFlt rtrevPi[20]; /* stationary frequencies for rtREV model */
1678 extern MrBFlt wagPi[20]; /* stationary frequencies for WAG model */
1679 extern MrBFlt cprevPi[20]; /* stationary frequencies for aacpREV model */
1680 extern MrBFlt vtPi[20]; /* stationary frequencies for VT model */
1681 extern MrBFlt blosPi[20]; /* stationary frequencies for Blosum62 model */
1682 extern MrBFlt lgPi[20]; /* stationary frequencies for LG model */
1684 #if defined (PRINT_DUMP)
1685 FILE *dumpFile; /* for debugging logs */
1688 #if defined (MPI_ENABLED)
1689 extern int proc_id; /* process ID (0, 1, ..., num_procs-1) */
1690 extern int num_procs; /* number of active processors */
1691 extern MrBFlt myStateInfo[7]; /* likelihood/prior/heat/ran/moveInfo vals of me */
1692 extern MrBFlt partnerStateInfo[7]; /* likelihood/prior/heat/ran/moveInfo vals of partner */
1695 #if defined (FAST_LOG)
1696 extern CLFlt scalerValue[];
1697 extern CLFlt logValue[];
1700 #endif /* __BAYES_H__ */