1 #ifndef __BAYES_H__ 2 #define __BAYES_H__ 3 4 #ifdef HAVE_CONFIG_H 5 # include "config.h" 6 # define VERSION_NUMBER PACKAGE_VERSION 7 #elif !defined (XCODE_VERSION) /* some defaults that would otherwise be guessed by configure */ 8 # define PACKAGE_NAME "MrBayes" 9 # define PACKAGE_VERSION "3.2.7" 10 # define HOST_CPU "x86_64" 11 # define VERSION_NUMBER PACKAGE_VERSION 12 # undef HAVE_LIBREADLINE 13 # define UNIX_VERSION 1 14 # define SSE_ENABLED 1 15 # undef AVX_ENABLED 16 # undef FMA_ENABLED 17 # undef MPI_ENABLED 18 # undef BEAGLE_ENABLED 19 #endif 20 21 #ifdef HAVE_UNISTD_H 22 #define _XOPEN_SOURCE 23 #include <unistd.h> 24 #endif 25 26 #include <assert.h> 27 #include <ctype.h> 28 #include <float.h> 29 #include <limits.h> 30 #include <math.h> 31 #include <memory.h> 32 #include <stdio.h> 33 #include <stdlib.h> 34 #include <string.h> 35 #include <stdarg.h> 36 #include <time.h> 37 38 /* Set SSE_ENABLED if SSE SIMD extensions available. */ 39 #ifdef HAVE_SSE 40 #define SSE_ENABLED 41 #endif 42 43 /* Set AVX_ENABLED if AVX SIMD extensions available. */ 44 #ifdef HAVE_AVX 45 #define AVX_ENABLED 46 #endif 47 48 /* Set FMA_ENABLED if FMA SIMD extensions available. */ 49 #if defined(HAVE_FMA3) || defined(HAVE_FMA4) 50 #define FMA_ENABLED 51 #endif 52 53 /* Set COMPLETIONMATCHES if we have the readline library */ 54 #ifdef HAVE_LIBREADLINE 55 #define COMPLETIONMATCHES 56 #endif 57 58 #if defined (MPI_ENABLED) 59 #include "mpi.h" 60 #endif 61 62 #if defined (BEAGLE_ENABLED) 63 #include "libhmsbeagle/beagle.h" 64 #endif 65 66 #if !defined (UNIX_VERSION) && !defined (WIN_VERSION) && !defined (MAC_VERSION) 67 # ifdef __MWERKS__ 68 # define MAC_VERSION 69 # elif defined __APPLE__ 70 # define MAC_VERSION 71 # elif defined __unix__ 72 # define UNIX_VERISON 73 # else 74 # define WIN_VERSION 75 # endif 76 #endif 77 78 #if defined (WIN_VERSION) 79 # include <windows.h> 80 # include <winbase.h> 81 #endif 82 83 /* Previous problems with bitfield operations may have been caused by several things. One potential 84 problem has been that MrBayes has used signed ints or signed longs for the bit operations, and 85 not all compilers perform bitfield operations on signed ints and longs in the expected manner 86 for unsigned ints and longs. Actually, the ANSI standard does not specify how to implement 87 bit operations on signed operands. 88 89 Another problem is that some bitshift operations in MrBayes have been performed on an integer 90 constant "1". This is probably interpreted by most compilers as a signed int and the result of 91 the bitshift operation is therefore also a signed int. When a signed int has had a different 92 number of bits than the type specified as SafeLong, this may also have resulted in errors on 93 some systems. 94 95 Both of these problems are fixed now by using unsigned longs for the bitfield operations and 96 by setting the "1" value needed by some bitfield functions using a variable defined as the same 97 type as the bitfield containers themselves. Furthermore, I have separated out the cases where 98 a signed long is required or is more appropriate than the container type of bitfields. 99 100 FR 2013-07-06 101 */ 102 typedef unsigned long BitsLong; 103 typedef long RandLong; 104 105 #define MRBFLT_MAX DBL_MAX /* maximum possible value that can be stored in MrBFlt */ 106 #define MRBFLT_MIN DBL_MIN /* minimum possible value that can be stored in MrBFlt */ 107 #define MRBFLT_NEG_MAX (-DBL_MAX) /* maximum possible negative value that can be stored in MrBFlt */ 108 typedef double MrBFlt; /* double used for parameter values and generally for floating point values, 109 if set to float MPI would not work becouse of use MPI_DOUBLE */ 110 typedef float CLFlt; /* single-precision float used for cond likes (CLFlt) to increase speed and reduce memory requirement */ 111 /* set CLFlt to double if you want increased precision */ 112 /* NOTE: CLFlt = double not compatible with SSE_ENABLED, AVX_ENABLED or FMA_ENABLED */ 113 114 /* 115 * Make sure we define SIMD instruction flags in a stepwise manner. That is, if we have FMA, make sure we have AVX; 116 * if we have AVX, make sure we have SSE. 117 */ 118 #if defined (FMA_ENABLED) 119 # if !defined (AVX_ENABLED) 120 # define AVX_ENABLED 121 # endif 122 #endif 123 124 #if defined (AVX_ENABLED) 125 # if !defined (SSE_ENABLED) 126 # define SSE_ENABLED 127 # endif 128 #endif 129 130 /* Likewise, if the user has asked for SIMD instructions to be disabled, do this 131 * too in a stepwise manner */ 132 #ifdef DISABLE_SSE 133 #undef SSE_ENABLED 134 #undef AVX_ENABLED 135 #undef FMA_ENABLED 136 #endif 137 138 #ifdef DISABLE_AVX 139 #undef AVX_ENABLED 140 #undef FMA_ENABLED 141 #endif 142 143 #ifdef DISABLE_FMA 144 #undef FMA_ENABLED 145 #endif 146 147 /* Define compiler for appropriate SIMD vector data alignment and free operations */ 148 #if defined (SSE_ENABLED) 149 # if defined (WIN_VERSION) 150 # define MS_VCPP_SIMD 151 # else 152 # define GCC_SIMD 153 # undef ICC_SIMD 154 # endif 155 #endif 156 157 /* Include header for SSE code */ 158 #if defined (SSE_ENABLED) 159 # include <xmmintrin.h> 160 #endif 161 162 /* Include header for AVX code */ 163 #if defined (AVX_ENABLED) 164 # include <immintrin.h> 165 #endif 166 167 #if defined (FMA_ENABLED) 168 # include <immintrin.h> 169 #endif 170 171 /* Define vector code constants */ 172 #define VEC_NONE 0 173 #define VEC_SSE 1 174 #define VEC_AVX 2 175 #define VEC_FMA 3 176 177 /* For comparing floating points: two values are the same if the absolute difference is less then 178 this value. 179 */ 180 #ifndef ETA 181 #define ETA (1E-30) 182 #endif 183 184 /* This defined DEBUG() is not used anywhere 185 #if defined (DEBUGOUTPUT) 186 #define DEBUG(fmt, arg) printf("%s:%d ",__FILE__,__LINE__); printf(fmt,arg); 187 #endif 188 */ 189 190 /* TEMPSTRSIZE determines size of temporary sprintf buffer (for SafeSprintf) */ 191 /* A patch was sent in by Allen Smith for SafeSprintf, but I could not get 192 it compiled on SGI IRIX 6.5 (too old?) with _xpg5_vsnprintf undefined. 193 The code below is a hack so SafeSprintf never has to reallocate memory. 194 */ 195 #ifdef __sgi 196 #define TEMPSTRSIZE 1000 197 #else 198 #define TEMPSTRSIZE 200 199 #endif 200 201 /* NO_ERROR is defined in bayes.h (as 0) and also in WinError.h (as 0L) 202 ERROR is defined in bayes.h (as 1) and also in WinGDI.h (as 0). we use the bayes.h value */ 203 #undef NO_ERROR 204 #undef ERROR 205 #define NO_ERROR 0 206 #define ERROR 1 207 #define NO_ERROR_QUIT 2 208 #define ABORT 3 209 #define SKIP_COMMAND 4 210 211 #undef FALSE 212 #undef TRUE 213 #define FALSE 0 214 #define TRUE 1 215 216 #define NO 0 217 #define YES 1 218 219 #define UP 0 220 #define DOWN 1 221 222 #define UPPER 0 223 #define MIDDLE 1 224 #define LOWER 2 225 226 #define NONINTERACTIVE 0 227 #define INTERACTIVE 1 228 229 #define STANDARD_USER 1 230 #define DEVELOPER 3 231 232 #define DIFFERENT 0 233 #define SAME 1 234 #define CONSISTENT_WITH 2 235 236 #define LINETERM_UNIX 0 237 #define LINETERM_MAC 1 238 #define LINETERM_DOS 2 239 240 #define SCREENWIDTH 60 241 #define SCREENWIDTH2 61 242 243 #define AVGSTDDEV 0 244 #define MAXSTDDEV 1 245 246 #define NONE 0 247 #define DNA 1 248 #define RNA 2 249 #define PROTEIN 3 250 #define RESTRICTION 4 251 #define STANDARD 5 252 #define MIXED 6 253 #define CONTINUOUS 7 254 255 #define AAMODEL_POISSON 0 256 #define AAMODEL_JONES 1 257 #define AAMODEL_DAY 2 258 #define AAMODEL_MTREV 3 259 #define AAMODEL_MTMAM 4 260 #define AAMODEL_WAG 5 261 #define AAMODEL_RTREV 6 262 #define AAMODEL_CPREV 7 263 #define AAMODEL_VT 8 264 #define AAMODEL_BLOSUM 9 265 #define AAMODEL_LG 10 266 #define AAMODEL_EQ 11 267 #define AAMODEL_GTR 12 /* aa models with free parameters must be listed last */ 268 269 #define NUCMODEL_4BY4 0 270 #define NUCMODEL_DOUBLET 1 271 #define NUCMODEL_CODON 2 272 #define NUCMODEL_AA 3 273 274 #define NST_MIXED -1 /* anything other than 1, 2, or 6 */ 275 276 #define MISSING 100000000 277 #define GAP 100000001 278 279 #define UNORD 0 280 #define ORD 1 281 #define DOLLO 2 282 #define IRREV 3 283 284 #define IN_CMD 0 285 #define IN_FILE 1 286 287 #define NOTHING 0 288 #define COMMAND 1 289 #define PARAMETER 2 290 #define EQUALSIGN 3 291 #define COLON 4 292 #define SEMICOLON 5 293 #define COMMA 6 294 #define POUNDSIGN 7 295 #define QUESTIONMARK 8 296 #define DASH 9 297 #define LEFTPAR 10 298 #define RIGHTPAR 11 299 #define LEFTCOMMENT 12 300 #define RIGHTCOMMENT 13 301 #define ALPHA 14 302 #define NUMBER 15 303 #define RETURNSYMBOL 16 304 #define ASTERISK 17 305 #define BACKSLASH 18 306 #define FORWARDSLASH 19 307 #define EXCLAMATIONMARK 20 308 #define PERCENT 21 309 #define QUOTATIONMARK 22 310 #define WEIRD 23 311 #define UNKNOWN_TOKEN_TYPE 24 312 #define LEFTCURL 25 313 #define RIGHTCURL 26 314 #define DOLLAR 27 315 #define AMPERSAND 28 316 #define VERTICALBAR 29 317 318 #define MAX_Q_RATE 100.0f 319 #define MIN_SHAPE_PARAM 0.00001f 320 #define MAX_SHAPE_PARAM 100.0f 321 #define MAX_SITE_RATE 10.0f 322 #define MAX_RATE_CATS 20 323 #define MAX_RATE_CATS_SQUARED 400 324 #define BRLENS_MIN 0.00000001f // 1E-8f 325 #define BRLENS_MAX 100.0f 326 /* BRLENS_MIN must be bigger than TIME_MIN */ 327 #define TIME_MIN 1.0E-11f 328 #define TIME_MAX 100.0f 329 #define RELBRLENS_MIN 0.00000001f // 1E-8f 330 #define RELBRLENS_MAX 100.0f 331 #define KAPPA_MIN 0.001f 332 #define KAPPA_MAX 1000.0f 333 #define GROWTH_MIN 0.000001f 334 #define GROWTH_MAX 1000000.0f 335 #define RATE_MIN 0.000001f 336 #define RATE_MAX 100.0f 337 #define CPPRATEMULTIPLIER_MIN 0.001f 338 #define CPPRATEMULTIPLIER_MAX 1000.0f 339 #define SYMPI_MIN 0.000001f 340 #define SYMPI_MAX 100.0f 341 #define ALPHA_MIN 0.0001f 342 #define ALPHA_MAX 10000.0f 343 #define DIR_MIN 0.000001f 344 #define PI_MIN 0.000001f 345 #define OFFSETEXPLAMBDA_MIN 0.000001f 346 #define OFFSETEXPLAMBDA_MAX 100000.0f 347 #define TREEHEIGHT_MIN 0.00000001f 348 #define TREEHEIGHT_MAX 1000.0f 349 #define TREEAGE_MIN 0.00000001f 350 #define TREEAGE_MAX 1000000.0f 351 #define CPPLAMBDA_MIN 0.00001f 352 #define CPPLAMBDA_MAX 100.0f 353 #define TK02VAR_MIN 0.000001f 354 #define TK02VAR_MAX 10000.0f 355 #define IGRVAR_MIN 0.000001f 356 #define IGRVAR_MAX 10000.0f 357 #define MIXEDVAR_MIN 0.000001f 358 #define MIXEDVAR_MAX 10000.0f 359 #define OMEGA_MIN 0.001f 360 #define OMEGA_MAX 1000.0f 361 362 #define POS_MIN 1E-25f 363 #define POS_MAX 1E25f 364 #define POS_INFINITY 1E25f 365 #define NEG_INFINITY -1000000.0f 366 #define POSREAL_MIN 1E-25f 367 #define POSREAL_MAX 1E25f 368 369 #define CMD_STRING_LENGTH 100000 370 371 #define pos(i,j,n) ((i)*(n)+(j)) 372 373 #define NUM_ALLOCS 100 374 375 #define ALLOC_MATRIX 0 376 #define ALLOC_CHARINFO 2 377 #define ALLOC_CHARSETS 3 378 #define ALLOC_TAXA 4 379 #define ALLOC_TMPSET 5 380 #define ALLOC_PARTITIONS 6 381 #define ALLOC_PARTITIONVARS 7 382 #define ALLOC_TAXASETS 8 383 #define ALLOC_CONSTRAINTS 9 384 #define ALLOC_USERTREE 10 385 #define ALLOC_SUMTPARAMS 11 386 #define ALLOC_TERMSTATE 12 387 #define ALLOC_ISPARTAMBIG 13 388 #define ALLOC_AVAILNODES 25 389 #define ALLOC_AVAILINDICES 26 390 #define ALLOC_CURLNL 28 391 #define ALLOC_CURLNPR 29 392 #define ALLOC_CHAINID 30 393 #define ALLOC_PARAMS 31 394 #define ALLOC_TREE 32 395 #define ALLOC_NODES 33 396 #define ALLOC_LOCTAXANAMES 34 397 #define ALLOC_COMPMATRIX 39 398 #define ALLOC_NUMSITESOFPAT 40 399 #define ALLOC_COMPCOLPOS 41 400 #define ALLOC_COMPCHARPOS 42 401 #define ALLOC_ORIGCHAR 43 402 #define ALLOC_PARAMVALUES 46 403 #define ALLOC_MCMCTREES 47 404 #define ALLOC_MOVES 48 405 #define ALLOC_PRELIKES 52 406 #define ALLOC_SITEJUMP 54 407 #define ALLOC_MARKOVTIS 55 408 #define ALLOC_RATEPROBS 56 409 #define ALLOC_STDTYPE 57 410 #define ALLOC_PACKEDTREES 58 411 #define ALLOC_SUMPSTRING 62 412 #define ALLOC_SUMPINFO 63 413 #define ALLOC_SWAPINFO 64 414 #define ALLOC_SYMPIINDEX 65 415 #define ALLOC_POSSELPROBS 66 416 #define ALLOC_PBF 68 417 #define ALLOC_LOCALTAXONCALIBRATION 69 418 #define ALLOC_SPR_PARSSETS 72 419 #define ALLOC_PFCOUNTERS 74 420 #define ALLOC_FILEPOINTERS 75 421 #define ALLOC_STATS 76 422 #define ALLOC_DIAGNTREE 77 423 #define ALLOC_USEDMOVES 82 424 #define ALLOC_MODEL 83 425 #define ALLOC_STDSTATEFREQS 84 426 #define ALLOC_PRINTPARAM 85 427 #define ALLOC_TREELIST 86 428 #define ALLOC_TFILEPOS 87 429 #define ALLOC_BEST 88 430 #define ALLOC_SPECIESPARTITIONS 89 431 #define ALLOC_SS 90 432 433 #define LINKED 0 434 #define UNLINKED 1 435 436 /*paramType*/ 437 #define NUM_LINKED 32 438 #define P_TRATIO 0 439 #define P_REVMAT 1 440 #define P_OMEGA 2 441 #define P_PI 3 442 #define P_SHAPE 4 443 #define P_PINVAR 5 444 #define P_CORREL 6 445 #define P_SWITCH 7 446 #define P_RATEMULT 8 447 #define P_TOPOLOGY 9 448 #define P_BRLENS 10 449 #define P_SPECRATE 11 450 #define P_EXTRATE 12 451 #define P_FOSLRATE 13 452 #define P_POPSIZE 14 453 #define P_AAMODEL 15 454 #define P_BRCORR 16 455 #define P_BRSIGMA 17 456 #define P_GROWTH 18 457 #define P_CPPMULTDEV 19 458 #define P_CPPRATE 20 459 #define P_CPPEVENTS 21 460 #define P_TK02VAR 22 461 #define P_TK02BRANCHRATES 23 462 #define P_IGRVAR 24 463 #define P_IGRBRANCHRATES 25 464 #define P_CLOCKRATE 26 465 #define P_SPECIESTREE 27 466 #define P_GENETREERATE 28 467 #define P_MIXEDVAR 29 468 #define P_MIXEDBRCHRATES 30 469 #define P_MIXTURE_RATES 31 470 /* NOTE: If you add another parameter, change NUM_LINKED */ 471 472 // #define CPPm 0 /* CPP rate multipliers */ 473 // #define CPPi 1 /* CPP independent rates */ 474 #define RCL_TK02 0 475 #define RCL_IGR 1 /* type of mixed relaxed clock model */ 476 477 #define MAX_NUM_USERTREES 200 /* maximum number of user trees MrBayes will read */ 478 #define MAX_CHAINS 256 /* maximum numbder of chains you can run actually only half of it becouse of m->lnLike[MAX_CHAINS] */ 479 480 // #define PARAM_NAME_SIZE 400 481 482 typedef void * VoidPtr; 483 typedef int (*CmdFxn)(void); 484 typedef int (*ParmFxn)(char *, char *); 485 486 /* typedef for a ln prior prob fxn */ 487 typedef MrBFlt (*LnPriorProbFxn)(MrBFlt val, MrBFlt *priorParams); 488 489 /* typedef for a ln prior prob ratio fxn */ 490 typedef MrBFlt (*LnPriorRatioFxn)(MrBFlt newVal, MrBFlt oldVal, MrBFlt *priorParams); 491 492 typedef struct 493 { 494 MrBFlt sum; /* sum of standard deviations */ 495 MrBFlt max; /* maximum standard deviation */ 496 MrBFlt numPartitions; 497 MrBFlt numSamples; 498 MrBFlt avgStdDev; 499 MrBFlt **pair; 500 } STATS; 501 502 /* enumeration for calibration prior type */ 503 enum CALPRIOR 504 { 505 unconstrained, 506 fixed, 507 uniform, 508 offsetExponential, 509 truncatedNormal, 510 logNormal, 511 offsetLogNormal, 512 standardGamma, 513 offsetGamma 514 }; 515 516 enum ConstraintType 517 { 518 PARTIAL, 519 NEGATIVE, 520 HARD 521 }; 522 523 enum CodingType 524 { 525 ALL = 0, 526 NOABSENCESITES = 1, 527 NOPRESENCESITES = 2, 528 VARIABLE = 3, 529 NOSINGLETONPRESENCE = 4, 530 NOSINGLETONABSENCE = 8, 531 NOSINGLETONS = 12, 532 INFORMATIVE = 15 533 }; 534 535 /* typedef for calibration */ 536 typedef struct calibration 537 { 538 char name[100]; 539 enum CALPRIOR prior; 540 MrBFlt priorParams[3]; 541 LnPriorProbFxn LnPriorProb; 542 LnPriorRatioFxn LnPriorRatio; 543 MrBFlt min; 544 MrBFlt max; 545 } 546 Calibration; 547 548 /* typedef for tree (topology) list element */ 549 typedef struct element 550 { 551 struct element *next; 552 int *order; 553 } TreeListElement; 554 555 /* typedef for list of trees (topologies) */ 556 typedef struct 557 { 558 TreeListElement *first; 559 TreeListElement *last; 560 } TreeList; 561 562 /* typedef for packed tree */ 563 typedef struct 564 { 565 int *order; 566 MrBFlt *brlens; 567 } PackedTree; 568 569 /* typedef for binary tree node */ 570 /* NOTE: Any variable added here must also be copied in CopyTrees */ 571 typedef struct node 572 { 573 char *label; /*!< name of node if tip */ 574 struct node *left, *right, *anc; /*!< pointers to adjacent nodes */ 575 int memoryIndex; /*!< memory index (do not change) */ 576 int index; /*!< index to node (0 to numLocalTaxa for tips) */ 577 int upDateCl; /*!< cond likes need update? */ 578 int upDateTi; /*!< transition probs need update? */ 579 int isLocked; /*!< is node locked? */ 580 int lockID; /*!< id of lock */ 581 int isDated; /*!< is node dated (calibrated)? */ 582 int marked, x, y; /*!< scratch variables */ 583 MrBFlt d; /*!< scratch variable */ 584 BitsLong *partition; /*!< pointer to bitfield describing splits */ 585 MrBFlt length; /*!< length of pending branch */ 586 MrBFlt nodeDepth; /*!< node depth (height) */ 587 MrBFlt age; /*!< age of node */ 588 Calibration *calibration; /*!< pointer to calibration data */ 589 } 590 TreeNode; 591 592 /* typedef for binary tree */ 593 typedef struct 594 { 595 char name[100]; /*!< name of tree */ 596 int memNodes; /*!< number of allocated nodes (do not exceed!) */ 597 int nNodes; /*!< number of nodes in tree (including lower root in rooted trees) */ 598 int nIntNodes; /*!< number of interior nodes in tree (excluding lower root in rooted trees) */ 599 int isRooted; /*!< is tree rooted? */ 600 int isClock; /*!< is tree clock? */ 601 int isCalibrated; /*!< is tree calibrated? */ 602 int nRelParts; /*!< number of relevant partitions */ 603 int *relParts; /*!< pointer to relevant partitions */ 604 int checkConstraints; /*!< does tree have constraints? */ 605 int nConstraints; /*!< number of constraints */ 606 int *constraints; /*!< pointer to constraints */ 607 int nLocks; /*!< number of constrained (locked) nodes */ 608 TreeNode **allDownPass; /*!< downpass array of all nodes */ 609 TreeNode **intDownPass; /*!< downpass array of interior nodes (including upper but excluding lower root in rooted trees) */ 610 #if defined (BEAGLE_V3_ENABLED) 611 int levelPassEnabled; /*!< are we also doing a level-order traversal? */ 612 TreeNode **intDownPassLevel; /*!< level order downpass array of interior nodes (including upper but excluding lower root in rooted trees) */ 613 #endif 614 TreeNode *root; /*!< pointer to root (lower root in rooted trees) */ 615 TreeNode *nodes; /*!< array containing the nodes */ 616 BitsLong *bitsets; /*!< pointer to bitsets describing splits */ 617 BitsLong *flags; /*!< pointer to cond like flags */ 618 int fromUserTree; /*!< YES is set for the trees whoes branch lengths are set from user tree(as start tree or fix branch length prior), NO otherwise */ 619 } 620 Tree; 621 622 /* typedef for node in polytomous tree */ 623 typedef struct pNode 624 { 625 char label[100]; /*!< name of node if terminal */ 626 struct pNode *left, *sib, *anc; /*!< pointers to adjacent nodes */ 627 int x, y, mark; /*!< scratch variables */ 628 int partitionIndex; /*!< partition index in sumt (scratch) */ 629 int index; /*!< index of node (if < numLocalTaxa = local taxon index) */ 630 int memoryIndex; /*!< immutable index of memory position */ 631 int isLocked; /*!< is the node locked? */ 632 int lockID; /*!< id of lock */ 633 int isDated; /*!< is node dated? */ 634 MrBFlt length; /*!< age of node */ 635 MrBFlt depth; /*!< depth (height) of node */ 636 MrBFlt age; /*!< age of node */ 637 MrBFlt support, f; /*!< scratch variables */ 638 BitsLong *partition; /*!< pointer to partition (split) bitset */ 639 Calibration *calibration; /*!< pointer to dating of node */ 640 } 641 PolyNode; 642 643 /* typedef for polytomous tree */ 644 typedef struct 645 { 646 char name[100]; /*!< name of tree */ 647 int memNodes; /*!< number of allocated nodes; do not exceed! */ 648 int nNodes; /*!< number of nodes in tree */ 649 int nIntNodes; /*!< number of interior nodes in tree */ 650 PolyNode **allDownPass; /*!< downpass array over all nodes */ 651 PolyNode **intDownPass; /*!< downpass array over interior nodes */ 652 PolyNode *root; /*!< pointer to root (lower for rooted trees */ 653 PolyNode *nodes; /*!< array holding the tree nodes */ 654 BitsLong *bitsets; /*!< bits describing partitions (splits) */ 655 int nBSets; /*!< number of effective branch length sets */ 656 int nESets; /*!< number of breakpoint rate sets */ 657 char **bSetName; /*!< names of effective branch length sets */ 658 char **eSetName; /*!< names of breakpoint rate sets */ 659 int **nEvents; /*!< number of branch events of bp rate set */ 660 MrBFlt ***position; /*!< position of branch events */ 661 MrBFlt ***rateMult; /*!< parameter of branch events */ 662 MrBFlt **effectiveBrLen; /*!< effective branch lengths of ebl set */ 663 int brlensDef; /*!< are brlens defined ? */ 664 int isRooted; /*!< is tree rooted? */ 665 int isClock; /*!< is tree clock? */ 666 int isCalibrated; /*!< is tree calibrated? */ 667 int isRelaxed; /*!< is tree relaxed? */ 668 MrBFlt clockRate; /*!< clock rate */ 669 int popSizeSet; /*!< does tree have a population size set? */ 670 MrBFlt *popSize; /*!< the population size */ 671 char *popSizeSetName; /*!< name of the population size set */ 672 } 673 PolyTree; 674 675 /* struct for holding model parameter info for the mcmc run */ 676 typedef struct param 677 { 678 int index; /* index to the parameter (0, 1, 2, ...) */ 679 int paramType; /* the type of the parameter */ 680 int paramId; /* unique ID for parameter x prior combination */ 681 MrBFlt *values; /* main values of parameter */ 682 MrBFlt *subValues; /* subvalues of parameter */ 683 int *intValues; /* integer values (model index/growth fxn) */ 684 int nValues; /* number of values */ 685 int nSubValues; /* number of subvalues */ 686 int nIntValues; /* number of intvalues */ 687 MrBFlt min; /* minimum value of parameter */ 688 MrBFlt max; /* maximum value of parameter */ 689 int *relParts; /* pointer to relevant divisions */ 690 int nRelParts; /* number of relevant divisions */ 691 int upDate; /* update flag (for copying) */ 692 struct param **subParams; /* pointers to subparams (for topology) */ 693 int nSubParams; /* number of subparams */ 694 Tree **tree; /* pointer to tree ptrs (for brlens & topology) */ 695 int treeIndex; /* index to first tree in mcmcTree */ 696 int hasBinaryStd; /* has binary standard chars */ 697 int *sympiBsIndex; /* pointer to sympi bsIndex (std chars) */ 698 int *sympinStates; /* pointer to sympi nStates (std chars) */ 699 int *sympiCType; /* pointer to sympi cType (std chars) */ 700 int nSympi; /* number of sympis */ 701 int printParam; /* whether parameter should be printed */ 702 int nPrintSubParams; /* number of subparams that should be printed */ 703 char *paramHeader; /* a string holding header for param values */ 704 char *name; /* string holding name of parameter */ 705 char *paramTypeName; /* pointer to description of parameter type */ 706 int checkConstraints; /* is tree parameter constrained? */ 707 int fill; /* flags whether the parameter should be filled */ 708 int nStdStateFreqs; /* number of std state frequencies */ 709 MrBFlt *stdStateFreqs; /* pointer to std state frequencies */ 710 int **nEvents; /* number of branch events for Cpp model */ 711 /* nEvents[0..2*numCains][0..numNodes=2*numTaxa] */ 712 MrBFlt ***position; /* event positions for Cpp relaxed clock model */ 713 MrBFlt ***rateMult; /* rate multipliers for Cpp relaxed clock model */ 714 int affectsLikelihood; /* does parameter directly influence likelihood? */ 715 MrBFlt* priorParams; /* pointer to the prior parameters */ 716 LnPriorProbFxn LnPriorProb; /* ln prior prob function */ 717 LnPriorRatioFxn LnPriorRatio; /* ln prior prob ratio function */ 718 } Param; 719 720 /* parameter ID values */ 721 /* identifies unique model parameter x prior combinations */ 722 #define TRATIO_DIR 1 723 #define TRATIO_FIX 2 724 #define REVMAT_DIR 3 725 #define REVMAT_FIX 4 726 #define OMEGA_DIR 5 727 #define OMEGA_FIX 6 728 #define SYMPI_UNI 7 729 #define SYMPI_UNI_MS 8 730 #define SYMPI_EXP 9 731 #define SYMPI_EXP_MS 10 732 #define SYMPI_FIX 11 733 #define SYMPI_FIX_MS 12 734 #define SYMPI_EQUAL 13 735 #define PI_DIR 14 736 #define PI_USER 15 737 #define PI_EMPIRICAL 16 738 #define PI_EQUAL 17 739 #define PI_FIXED 18 740 #define SHAPE_UNI 19 741 #define SHAPE_EXP 20 742 #define SHAPE_FIX 21 743 #define PINVAR_UNI 22 744 #define PINVAR_FIX 23 745 #define CORREL_UNI 24 746 #define CORREL_FIX 25 747 #define SWITCH_UNI 26 748 #define SWITCH_EXP 27 749 #define SWITCH_FIX 28 750 #define RATEMULT_DIR 29 751 #define RATEMULT_FIX 30 752 #define TOPOLOGY_NCL_UNIFORM 31 753 #define TOPOLOGY_NCL_CONSTRAINED 32 754 #define TOPOLOGY_NCL_FIXED 33 755 #define TOPOLOGY_NCL_UNIFORM_HOMO 34 756 #define TOPOLOGY_NCL_UNIFORM_HETERO 35 757 #define TOPOLOGY_NCL_CONSTRAINED_HOMO 36 758 #define TOPOLOGY_NCL_CONSTRAINED_HETERO 37 759 #define TOPOLOGY_NCL_FIXED_HOMO 38 760 #define TOPOLOGY_NCL_FIXED_HETERO 39 761 #define TOPOLOGY_CL_UNIFORM 40 762 #define TOPOLOGY_CL_CONSTRAINED 41 763 #define TOPOLOGY_CL_FIXED 42 764 #define TOPOLOGY_CCL_UNIFORM 43 765 #define TOPOLOGY_CCL_CONSTRAINED 44 766 #define TOPOLOGY_CCL_FIXED 45 767 #define TOPOLOGY_PARSIMONY_UNIFORM 46 768 #define TOPOLOGY_PARSIMONY_CONSTRAINED 47 769 #define TOPOLOGY_PARSIMONY_FIXED 48 770 #define BRLENS_UNI 49 771 #define BRLENS_EXP 50 772 #define BRLENS_GamDir 51 773 #define BRLENS_iGmDir 52 774 #define BRLENS_twoExp 53 775 #define BRLENS_FIXED 54 776 #define BRLENS_CLOCK_UNI 55 777 #define BRLENS_CLOCK_COAL 56 778 #define BRLENS_CLOCK_BD 57 779 #define BRLENS_CLOCK_FIXED 58 780 #define BRLENS_CLOCK_SPCOAL 59 781 #define BRLENS_CLOCK_FOSSIL 60 782 #define BRLENS_PARSIMONY 61 783 #define SPECRATE_UNI 62 784 #define SPECRATE_EXP 63 785 #define SPECRATE_FIX 64 786 #define EXTRATE_BETA 65 787 #define EXTRATE_FIX 66 788 #define FOSLRATE_BETA 67 789 #define FOSLRATE_FIX 68 790 #define POPSIZE_UNI 69 791 #define POPSIZE_GAMMA 70 792 #define POPSIZE_FIX 71 793 #define POPSIZE_NORMAL 72 794 #define POPSIZE_LOGNORMAL 73 795 #define AAMODEL_FIX 74 796 #define AAMODEL_MIX 75 797 #define GROWTH_UNI 76 798 #define GROWTH_EXP 77 799 #define GROWTH_FIX 78 800 #define GROWTH_NORMAL 79 801 #define OMEGA_BUD 80 802 #define OMEGA_BUF 81 803 #define OMEGA_BED 82 804 #define OMEGA_BEF 83 805 #define OMEGA_BFD 84 806 #define OMEGA_BFF 85 807 #define OMEGA_FUD 86 808 #define OMEGA_FUF 87 809 #define OMEGA_FED 88 810 #define OMEGA_FEF 89 811 #define OMEGA_FFD 90 812 #define OMEGA_FFF 91 813 #define OMEGA_ED 92 814 #define OMEGA_EF 93 815 #define OMEGA_FD 94 816 #define OMEGA_FF 95 817 #define OMEGA_10UUB 96 818 #define OMEGA_10UUF 97 819 #define OMEGA_10UEB 98 820 #define OMEGA_10UEF 99 821 #define OMEGA_10UFB 100 822 #define OMEGA_10UFF 101 823 #define OMEGA_10EUB 102 824 #define OMEGA_10EUF 103 825 #define OMEGA_10EEB 104 826 #define OMEGA_10EEF 105 827 #define OMEGA_10EFB 106 828 #define OMEGA_10EFF 107 829 #define OMEGA_10FUB 108 830 #define OMEGA_10FUF 109 831 #define OMEGA_10FEB 110 832 #define OMEGA_10FEF 111 833 #define OMEGA_10FFB 112 834 #define OMEGA_10FFF 113 835 #define CPPRATE_FIX 114 836 #define CPPRATE_EXP 115 837 #define CPPMULTDEV_FIX 116 838 #define TK02VAR_FIX 117 839 #define TK02VAR_EXP 118 840 #define TK02VAR_UNI 119 841 #define TOPOLOGY_RCL_UNIFORM 120 842 #define TOPOLOGY_RCL_CONSTRAINED 121 843 #define TOPOLOGY_RCL_FIXED 122 844 #define TOPOLOGY_RCCL_UNIFORM 123 845 #define TOPOLOGY_RCCL_CONSTRAINED 124 846 #define TOPOLOGY_RCCL_FIXED 125 847 #define TOPOLOGY_SPECIESTREE 126 848 #define CPPEVENTS 127 849 #define TK02BRANCHRATES 128 850 #define TOPOLOGY_FIXED 129 851 #define IGRVAR_FIX 130 852 #define IGRVAR_EXP 131 853 #define IGRVAR_UNI 132 854 #define IGRBRANCHRATES 133 855 #define CLOCKRATE_FIX 134 856 #define CLOCKRATE_NORMAL 135 857 #define CLOCKRATE_LOGNORMAL 136 858 #define CLOCKRATE_GAMMA 137 859 #define CLOCKRATE_EXP 138 860 #define SPECIESTREE_UNIFORM 139 861 #define GENETREERATEMULT_DIR 140 862 #define GENETREERATEMULT_FIX 141 863 #define REVMAT_MIX 142 864 #define MIXEDVAR_FIX 143 865 #define MIXEDVAR_EXP 144 866 #define MIXEDVAR_UNI 145 867 #define MIXEDBRCHRATES 146 868 #define MIXTURE_RATES 147 869 870 #if defined (BEAGLE_ENABLED) 871 #define MB_BEAGLE_SCALE_ALWAYS 0 872 #define MB_BEAGLE_SCALE_DYNAMIC 1 873 #if defined (_DEBUG) 874 #define MB_PRINT_DYNAMIC_RESCALE_FAIL_STAT 875 #endif 876 #endif 877 878 /* typedef for a MoveFxn */ 879 typedef int (MoveFxn)(Param *param, int chain, RandLong *seed, MrBFlt *lnPriorRatio, MrBFlt *lnProposalRatio, MrBFlt *mvp); 880 881 /* typedef for an ApplicFxn */ 882 typedef int (ApplicFxn)(Param *param); 883 884 /* typedef for an AutotuneFxn */ 885 typedef void (AutotuneFxn)(MrBFlt acceptanceRate, MrBFlt targetRate, int batch, MrBFlt *tuningParameter, MrBFlt minTuning, MrBFlt maxTuning); 886 887 /* struct holding info on each move type that the program handles */ 888 typedef struct 889 { 890 MoveFxn *moveFxn; /* pointer to the move function */ 891 ApplicFxn *isApplicable; /* pointer to function determining whether move is applicable to a parameter */ 892 int nApplicable; /* number of relevant params */ 893 int applicableTo[40]; /* pointer to ID of relevant params */ 894 char *name; /* name of the move type */ 895 char *shortName; /* abbreviated name of the move type */ 896 char *paramName; /* name of subparameter if complex parameter */ 897 int subParams; /* are we changing subparams (brlens of topol.) */ 898 char *tuningName[5]; /* name of tuning params */ 899 char *shortTuningName[5];/* short name of tuning params */ 900 MrBFlt relProposalProb; /* default relative proposal probability */ 901 int numTuningParams; /* number of tuning parameters */ 902 MrBFlt tuningParam[5]; /* default tuning parameters for the proposal */ 903 MrBFlt minimum[5]; /* minimum values for tuning params */ 904 MrBFlt maximum[5]; /* maximum values for tuning params */ 905 int parsimonyBased; /* this move is based on parsimony (YES/NO) */ 906 int level; /* user level of this move */ 907 AutotuneFxn *Autotune; /* pointer to the autotuning function */ 908 MrBFlt targetRate; /* default target acceptance rate for autotuning */ 909 } MoveType; 910 911 /* max number of move types */ 912 #define NUM_MOVE_TYPES 101 913 914 /* struct holding info on each move */ 915 /* Note: This allows several different moves to affect the same parameter */ 916 /* It also allows the same move to affect different parameters as before */ 917 /* This is also a good place to keep the proposal probs */ 918 typedef struct 919 { 920 char *name; /* name of the move */ 921 MoveType *moveType; /* pointer to the move type */ 922 MoveFxn *moveFxn; /* pointer to the move function */ 923 Param *parm; /* ptr to parameter the move applies to */ 924 MrBFlt *relProposalProb; /* the relative proposal probability */ 925 MrBFlt *cumProposalProb; /* the cumulative proposal probability */ 926 int *nAccepted; /* number of accepted moves */ 927 int *nTried; /* number of tried moves */ 928 int *nBatches; /* counter for autotuning rounds */ 929 int *nTotAccepted; /* total number of accepted moves */ 930 int *nTotTried; /* total number of tried moves */ 931 MrBFlt *targetRate; /* target acceptance rate for autotuning */ 932 MrBFlt *lastAcceptanceRate;/* acceptance rate in last complete batch */ 933 MrBFlt **tuningParam; /* tuning parameters for the move */ 934 } MCMCMove; 935 936 typedef int (*LikeDownFxn)(TreeNode *, int, int); 937 typedef int (*LikeRootFxn)(TreeNode *, int, int); 938 typedef int (*LikeScalerFxn)(TreeNode *, int, int); 939 typedef int (*LikeFxn)(TreeNode *, int, int, MrBFlt *, int); 940 typedef int (*TiProbFxn)(TreeNode *, int, int); 941 typedef int (*LikeUpFxn)(TreeNode *, int, int); 942 typedef int (*PrintAncStFxn)(TreeNode *, int, int); 943 typedef int (*StateCodeFxn) (int); 944 typedef int (*PrintSiteRateFxn) (TreeNode *, int, int); 945 typedef int (*PosSelProbsFxn) (TreeNode *, int, int); 946 typedef int (*SiteOmegasFxn) (TreeNode *, int, int); 947 948 typedef struct cmdtyp 949 { 950 int cmdNumber; 951 char *string; 952 int specialCmd; 953 CmdFxn cmdFxnPtr; 954 short numParms; 955 short parmList[50]; 956 int expect; 957 char *cmdDescription; 958 int cmdUse; 959 int hiding; 960 } 961 CmdType; 962 963 typedef struct parm 964 { 965 char *string; /* parameter name */ 966 char *valueList; /* list of values that could be input */ 967 ParmFxn fp; /* function pointer */ 968 } 969 ParmInfo, *ParmInfoPtr; 970 971 typedef struct model 972 { 973 int dataType; /* data type for partition */ 974 int nStates; /* number of states for this type of data */ 975 int codon[64]; /* gives protein ID for each codon */ 976 int codonNucs[64][3]; /* gives triplet for each codon */ 977 int codonAAs[64]; /* gives protein ID for implemented code */ 978 979 char nucModel[100]; /* nucleotide model used */ 980 char nst[100]; /* number of substitution types */ 981 char parsModel[100]; /* use the (so-called) parsimony model */ 982 char geneticCode[100]; /* genetic code used */ 983 int coding; /* type of patterns encoded */ 984 char codingString[100]; /* string describing type of patterns encoded */ 985 char ploidy[100]; /* ploidy level */ 986 char omegaVar[100]; /* type of omega variation model */ 987 char ratesModel[100]; /* rates across sites model */ 988 int numGammaCats; /* number of categories for gamma approximation */ 989 int numLnormCats; /* number of categories for lnorm approximation */ 990 int numMixtCats; /* number of components of rate mixture */ 991 char useGibbs[100]; /* flags whether Gibbs sampling of discrete gamma is used */ 992 int gibbsFreq; /* frequency of Gibbs resampling of discrete gamma */ 993 994 int numBetaCats; /* number of categories for beta approximation */ 995 int numM10GammaCats; /* number of cats for gamma approx (M10 model) */ 996 int numM10BetaCats; /* number of cats for beta approx (M10 model) */ 997 char covarionModel[100];/* use covarion model? (yes/no) */ 998 char augmentData[100]; /* should data be augmented */ 999 1000 char tRatioPr[100]; /* prior for ti/tv rate ratio */ 1001 MrBFlt tRatioFix; 1002 MrBFlt tRatioDir[2]; 1003 char revMatPr[100]; /* prior for GTR model */ 1004 MrBFlt revMatFix[6]; 1005 MrBFlt revMatDir[6]; 1006 MrBFlt revMatSymDir; /* prior for mixed GTR subspace model */ 1007 char aaModelPr[100]; /* prior for amino acid model */ 1008 char aaModel[100]; 1009 MrBFlt aaModelPrProbs[10]; 1010 char aaRevMatPr[100]; /* prior for aa GTR model */ 1011 MrBFlt aaRevMatFix[190]; 1012 MrBFlt aaRevMatDir[190]; 1013 char omegaPr[100]; /* prior for omega */ 1014 MrBFlt omegaFix; 1015 MrBFlt omegaDir[2]; 1016 char ny98omega1pr[100]; /* prior for class 1 omega (Ny98 model) */ 1017 MrBFlt ny98omega1Fixed; 1018 MrBFlt ny98omega1Beta[2]; 1019 char ny98omega3pr[100]; /* prior for class 3 omega (Ny98 model) */ 1020 MrBFlt ny98omega3Fixed; 1021 MrBFlt ny98omega3Uni[2]; 1022 MrBFlt ny98omega3Exp; 1023 char m3omegapr[100]; /* prior for all three omegas (M3 model) */ 1024 MrBFlt m3omegaFixed[3]; 1025 char m10betapr[100]; /* prior for omega variation (M10 model) */ 1026 char m10gammapr[100]; 1027 MrBFlt m10betaExp; 1028 MrBFlt m10betaUni[2]; 1029 MrBFlt m10betaFix[2]; 1030 MrBFlt m10gammaExp; 1031 MrBFlt m10gammaUni[2]; 1032 MrBFlt m10gammaFix[2]; 1033 char codonCatFreqPr[100]; /* prior for selection cat frequencies */ 1034 MrBFlt codonCatFreqFix[3]; 1035 MrBFlt codonCatDir[3]; 1036 char stateFreqPr[100]; /* prior for character state frequencies */ 1037 MrBFlt stateFreqsFix[200]; 1038 MrBFlt stateFreqsDir[200]; 1039 char stateFreqsFixType[100]; 1040 int numDirParams; 1041 char shapePr[100]; /* prior for gamma/lnorm shape parameter */ 1042 MrBFlt shapeFix; 1043 MrBFlt shapeUni[2]; 1044 MrBFlt shapeExp; 1045 char pInvarPr[100]; /* prior for proportion of invariable sites */ 1046 MrBFlt pInvarFix; 1047 MrBFlt pInvarUni[2]; 1048 char adGammaCorPr[100]; /* prior for correlation param of adGamma model */ 1049 MrBFlt corrFix; 1050 MrBFlt corrUni[2]; 1051 char covSwitchPr[100]; /* prior for switching rates of covarion model */ 1052 MrBFlt covswitchFix[2]; 1053 MrBFlt covswitchUni[2]; 1054 MrBFlt covswitchExp; 1055 char symPiPr[100]; /* prior for pi when unidentifiable states used */ 1056 MrBFlt symBetaFix; 1057 MrBFlt symBetaUni[2]; 1058 MrBFlt symBetaExp; 1059 char ratePr[100]; /* prior on rate for a partition */ 1060 MrBFlt ratePrDir; 1061 char generatePr[100]; /* prior on rate for a gene (one or more partitions) */ 1062 MrBFlt generatePrDir; 1063 char brownCorPr[100]; /* prior for correlation of Brownian model */ 1064 MrBFlt brownCorrFix; 1065 MrBFlt brownCorrUni[2]; 1066 char brownScalesPr[100]; /* prior for scales of Brownian model */ 1067 MrBFlt brownScalesFix; 1068 MrBFlt brownScalesUni[2]; 1069 MrBFlt brownScalesGamma[2]; 1070 MrBFlt brownScalesGammaMean; 1071 1072 char topologyPr[100]; /* prior for tree topology */ 1073 int topologyFix; /* user tree index for fixed topology */ 1074 int *activeConstraints; /* which constraints are active? */ 1075 int numActiveConstraints; 1076 int numActiveLocks; 1077 char brlensPr[100]; /* prior on branch lengths */ 1078 int brlensFix; /* user tree index for fixed brlens */ 1079 MrBFlt brlensUni[2]; 1080 MrBFlt brlensExp; 1081 MrBFlt brlens2Exp[2]; 1082 MrBFlt brlensDir[4]; 1083 MrBFlt brlensGamma[2]; 1084 char speciesTreeBrlensPr[100]; /* prior on branch lengths of species tree */ 1085 char unconstrainedPr[100]; /* prior on branch lengths if unconstrained */ 1086 char clockPr[100]; /* prior on branch if clock enforced */ 1087 char clockVarPr[100]; /* prior on clock rate variation (strict, cpp, tk02, igr, ...) */ 1088 char nodeAgePr[100]; /* prior on node depths (unconstrained, constraints) */ 1089 char speciationPr[100]; /* prior on speciation rate (net diversification) */ 1090 MrBFlt speciationFix; 1091 MrBFlt speciationUni[2]; 1092 MrBFlt speciationExp; 1093 char extinctionPr[100]; /* prior on relative extinction rate (turnover) */ 1094 MrBFlt extinctionFix; 1095 MrBFlt extinctionBeta[2]; 1096 char fossilizationPr[100]; /* prior on fossilization rate (sampling proportion) */ 1097 MrBFlt fossilizationFix; 1098 MrBFlt fossilizationBeta[2]; 1099 char sampleStrat[100]; /* taxon sampling strategy (for b-d process) */ 1100 int birthRateShiftNum; /* number of birth rate shifts */ 1101 MrBFlt birthRateShiftTime[100]; /* birth rate shifting times */ 1102 int deathRateShiftNum; /* number of death rate shifts */ 1103 MrBFlt deathRateShiftTime[100]; /* death rate shifting times */ 1104 int fossilSamplingNum; /* number of fossil sampling rate shiftings / slice sampling events */ 1105 MrBFlt fossilSamplingTime[100]; /* fossil sampling rate shifting times */ 1106 // MrBFlt *fossilSamplingProb; /* fossil slice sampling probs (rho_i) */ 1107 MrBFlt sampleProb; /* extant taxon sampling fraction (rho) */ 1108 Calibration treeAgePr; /* prior on tree age for uniform clock trees */ 1109 char clockRatePr[100]; /* prior on base substitution rate of tree for clock trees */ 1110 MrBFlt clockRateNormal[2]; 1111 MrBFlt clockRateLognormal[2]; 1112 MrBFlt clockRateGamma[2]; 1113 MrBFlt clockRateExp; 1114 MrBFlt clockRateFix; 1115 char popSizePr[100]; /* prior on population size */ 1116 MrBFlt popSizeFix; 1117 MrBFlt popSizeUni[2]; 1118 MrBFlt popSizeLognormal[2]; 1119 MrBFlt popSizeGamma[2]; 1120 MrBFlt popSizeNormal[2]; 1121 char popVarPr[100]; /* prior on pop. size variation across tree */ 1122 char growthPr[100]; /* prior on coalescence growth rate */ 1123 MrBFlt growthFix; 1124 MrBFlt growthUni[2]; 1125 MrBFlt growthExp; 1126 MrBFlt growthNorm[2]; 1127 char cppRatePr[100]; /* prior on CPP rate */ 1128 MrBFlt cppRateFix; 1129 MrBFlt cppRateExp; 1130 char cppMultDevPr[100]; /* prior on CPP rate multiplier Lognormal variance */ 1131 MrBFlt cppMultDevFix; 1132 char tk02varPr[100]; /* prior on TK02 lognormal rate variance */ 1133 MrBFlt tk02varFix; 1134 MrBFlt tk02varUni[2]; 1135 MrBFlt tk02varExp; 1136 char igrvarPr[100]; /* prior on IGR gamma distribution variance */ 1137 MrBFlt igrvarFix; 1138 MrBFlt igrvarUni[2]; 1139 MrBFlt igrvarExp; 1140 char mixedvarPr[100]; /* prior on mixed relaxed clock rate variance */ 1141 MrBFlt mixedvarFix; 1142 MrBFlt mixedvarUni[2]; 1143 MrBFlt mixedvarExp; 1144 1145 char tratioFormat[30]; /* format used to report tratio */ 1146 char revmatFormat[30]; /* format used to report revmat */ 1147 char ratemultFormat[30]; /* format used to report ratemult */ 1148 char treeFormat[30]; /* format used to report trees/topologies */ 1149 char inferAncStates[5]; /* should ancestral states be inferred (Yes/No)? */ 1150 char inferSiteOmegas[5]; /* should site omega vals be inferred (Yes/No)? */ 1151 char inferSiteRates[5]; /* should site rates be inferred (Yes/No)? */ 1152 char inferPosSel[5]; /* should site selection be inferred (Yes/No)? */ 1153 } Model, ModelParams; 1154 1155 typedef struct chain 1156 { 1157 int numGen; /* number of MCMC cycles */ 1158 int sampleFreq; /* frequency to sample chain */ 1159 int printFreq; /* frequency to print chain */ 1160 int swapFreq; /* frequency to attempt swap of states */ 1161 int numRuns; /* number of runs */ 1162 int numChains; /* number of chains */ 1163 int isSS; /* do we do Steppingstone Sampling */ 1164 int startFromPriorSS; /* If Yes SS is moving from Prior to Posterior */ 1165 int numStepsSS; /* Number of steps in SS */ 1166 int burninSS; /* Fixed burnin for SS */ 1167 MrBFlt alphaSS; /* Beta values are distributed according to quantiles of Beta(alphaSS,1.0) distribution */ 1168 int backupCheckSS; /* Frequency of checkpoints backup */ 1169 MrBFlt chainTemp; /* chain temperature */ 1170 int userDefinedTemps; /* should we use the users temperatures? */ 1171 MrBFlt userTemps[MAX_CHAINS]; /* user-defined chain temperatures */ 1172 char chainFileName[100]; /* chain file name for output */ 1173 int chainBurnIn; /* chain burn in length */ 1174 int numStartPerts; /* number of perturbations to starting tree */ 1175 char startTree[100]; /* starting tree for chain (current/random) */ 1176 char startParams[100]; /* starting values for chain (current/reset) */ 1177 int saveBrlens; /* should branch lengths be saved */ 1178 MrBFlt weightScheme[3]; /* percent chars to increase/decrease in weight */ 1179 int calcPbf; /* should we calculate the pseudo Bayes factor */ 1180 int pbfInitBurnin; /* initial burnin when calculating pseudo BF */ 1181 int pbfSampleFreq; /* sample frequency for pseudo BF */ 1182 int pbfSampleTime; /* how many cycles to calcualate site prob. */ 1183 int pbfSampleBurnin; /* burnin period for each site for pseudo BF */ 1184 int swapAdjacentOnly; /* whether we only swap adjacent temperatures */ 1185 int redirect; /* should output go to stdout */ 1186 int allChains; /* should stats be output for all chains? */ 1187 int numSwaps; /* number of swaps to try each time */ 1188 int mcmcDiagn; /* should mcmc diagnostics be output? */ 1189 int diagnFreq; /* mcmc diagnostics frequency */ 1190 int diagnStat; /* statistic to use for mcmc diagnostics */ 1191 int relativeBurnin; /* should a relative burnin be used ? */ 1192 MrBFlt burninFraction; /* the sample fraction to discard as burnin */ 1193 int allComps; /* top conv diagnosis for all pairs? */ 1194 MrBFlt minPartFreq; /* minimum partition frequency for conv diagn */ 1195 MrBFlt stopVal; /* top conv diagn value to reach before stopping */ 1196 int stopRule; /* use stop rule? */ 1197 STATS *stat; /* ptr to structs with mcmc diagnostics info */ 1198 Tree *dtree; /* pointing to tree used for conv diagnostics */ 1199 TreeList *treeList; /* vector of tree lists for saving trees */ 1200 int saveTrees; /* save tree samples for later removal? */ 1201 int stopTreeGen; /* generation after which no trees need be saved */ 1202 fpos_t *tFilePos; /* position for reading trees for removal */ 1203 int printMax; /* maximum number of chains to print */ 1204 int printAll; /* whether to print all or only cold chains */ 1205 int checkPoint; /* should we use check-pointing? */ 1206 int checkFreq; /* check-pointing frequency */ 1207 int runWithData; /* should we run with or without data? */ 1208 int orderTaxa; /* order taxa before printing tree to file? */ 1209 int append; /* order taxa before printing tree to file? */ 1210 int autotune; /* autotune tuning parameters of proposals ? */ 1211 int tuneFreq; /* autotuning frequency */ 1212 } Chain; 1213 1214 typedef struct modelinfo 1215 { 1216 /* General model information */ 1217 int dataType; /* data type for partition */ 1218 int nucModelId; /* structure of nucleotide model */ 1219 int nst; /* # substitution types */ 1220 int aaModelId; /* amino acid model type */ 1221 int parsModelId; /* is parsimony model used YES/NO */ 1222 1223 /* Specific model information */ 1224 int numRateCats; /* number of rate cats (1 if inapplic.) */ 1225 int numBetaCats; /* number of beta cats (1 if inapplic.) */ 1226 int numOmegaCats; /* number of omega cats (1 if inapplic.) */ 1227 int numTiCats; /* number of cats needing different tis */ 1228 int numModelStates; /* number of states including hidden ones */ 1229 int numStates; /* number of observable discrete states */ 1230 1231 /* Model parameter pointers */ 1232 Param *tRatio; /* ptr to tRatio used in model */ 1233 Param *revMat; /* ptr to revMat used in model */ 1234 Param *omega; /* ptr to omega used in model */ 1235 Param *stateFreq; /* ptr to statFreq used in model */ 1236 Param *mixtureRates; /* ptr to site rate mixture used in model */ 1237 Param *shape; /* ptr to shape used in model */ 1238 Param *pInvar; /* ptr to pInvar used in model */ 1239 Param *correlation; /* ptr to correlation used in model */ 1240 Param *switchRates; /* ptr to switchRates (off<->on) */ 1241 Param *rateMult; /* ptr to parition rateMult used in model */ 1242 Param *geneTreeRateMult; /* ptr to gene tree rateMult used in model */ 1243 Param *speciationRates; /* ptr to speciationRates used in model */ 1244 Param *extinctionRates; /* ptr to extinctionRates used in model */ 1245 Param *fossilizationRates; /* ptr to fossilizationRates */ 1246 Param *popSize; /* ptr to population size used in model */ 1247 Param *growthRate; /* ptr to growth rate used in model */ 1248 Param *topology; /* ptr to topology used in model */ 1249 Param *brlens; /* ptr to brlens (and tree) used in model */ 1250 Param *speciesTree; /* ptr to species tree used in model */ 1251 Param *aaModel; /* ptr to amino acid matrix used */ 1252 Param *cppMultDev; /* ptr to cpp ratemult lognormal variance */ 1253 Param *cppRate; /* ptr to CPP rate used in model */ 1254 Param *cppEvents; /* ptr to CPP events */ 1255 Param *tk02var; /* ptr to variance for TK02 relaxed clock */ 1256 Param *tk02BranchRates; /* ptr to branch rates for TK02 relaxed clock */ 1257 Param *igrvar; /* ptr to gamma var for IGR relaxed clock */ 1258 Param *igrBranchRates; /* ptr to branch rates for IGR relaxed clock*/ 1259 Param *mixedvar; /* ptr to var for mixed relaxed clock */ 1260 Param *mixedBrchRates; /* ptr to branch rates for mixed relaxed clock */ 1261 Param *clockRate; /* ptr to clock rate parameter */ 1262 1263 /* Information about characters and transformations */ 1264 int numChars; /* number of compressed characters */ 1265 int numUncompressedChars; /* number of uncompressed characters */ 1266 int numDummyChars; /* number of dummy characters */ 1267 int compMatrixStart; /* start column in compressed matrix */ 1268 int compMatrixStop; /* stop column in compressed matrix */ 1269 int compCharStart; /* start char among compressed chars */ 1270 int compCharStop; /* stop char among compressed chars */ 1271 int parsMatrixStart; /* start column in parsimony matrix */ 1272 int parsMatrixStop; /* stop collumn in parsimony matrix */ 1273 int nParsIntsPerSite; /* # parsimony ints per character */ 1274 int nCharsPerSite; /* number chars per site (eg 3 for codon) */ 1275 int rateProbStart; /* start of rate probs (for adgamma) */ 1276 1277 /* Variables for eigen decompositions */ 1278 int cijkLength; /* stores length of cijk vector */ 1279 int nCijkParts; /* stores number of cijk partitions (1 except for omega/covarion models) */ 1280 int upDateCijk; /* whether cijk vector needs to be updated */ 1281 1282 /* Variables for standard model */ 1283 int *tiIndex; /* index to trans probs for each compressed char*/ 1284 int *bsIndex; /* index to stat freqs for each compressed char */ 1285 int *nStates; /* # states of each compressed char */ 1286 int *cType; /* whether char is ord, unord or irrev */ 1287 int *weight; /* prior weight of each compressed char */ 1288 int isTiNeeded[20]; /* marks whether a trans prob matrix is needed */ 1289 1290 /* Gibbs sampling of gamma site rate parameters */ 1291 CLFlt ***catLike; /* likelihood for Gibbs sampling of gamma */ 1292 CLFlt ***catLnScaler; /* scaler for Gibbs sampling of gamma */ 1293 int gibbsGamma; /* flags whether Gibbs sampling of discrete gamma is used */ 1294 int gibbsFreq; /* frequency of Gibbs resampling of discrete gamma */ 1295 1296 /* Variables for parsimony sets and parsimony calculations */ 1297 MrBFlt parsTreeLength[MAX_CHAINS*2]; /* parsimony tree lengths for chains */ 1298 BitsLong **parsSets; /* parsimony sets */ 1299 int numParsSets; /* number of parsimony sets */ 1300 CLFlt *parsNodeLens; /* parsimony node lengths */ 1301 int numParsNodeLens; /* number of parsimony node lengths */ 1302 1303 /* Miscellaneous parameters */ 1304 int mark; /* scratch parameter */ 1305 int parsimonyBasedMove; /* is parsimony-based move used (YES/NO) */ 1306 1307 /* Variables for conditional likelihoods */ 1308 int upDateCl; /* flags whether update of cond likes needed */ 1309 int upDateAll; /* flags whether update of entire tree is needed*/ 1310 int *isPartAmbig; /* is tip partially ambiguous? */ 1311 int **termState; /* index arrays for terminal branch shortcuts */ 1312 CLFlt *invCondLikes; /* space for the invariable cond likes */ 1313 CLFlt **condLikes; /* space for the cond likes */ 1314 CLFlt **tiProbs; /* space for the ti probs */ 1315 CLFlt **scalers; /* space for the node and site scalers */ 1316 CLFlt **clP; /* handy pointers to cond likes for ti cats */ 1317 #if defined (SSE_ENABLED) 1318 __m128 **clP_SSE; /* handy pointers to cond likes, SSE version */ 1319 int numVecChars; /* number of compact SIMD vectors */ 1320 int numFloatsPerVec; /* number of floats per SIMD vector */ 1321 CLFlt *lnL_Vec; /* temp storage for log site likes */ 1322 CLFlt *lnLI_Vec; /* temp storage for log site invariable likes */ 1323 #if defined (AVX_ENABLED) 1324 __m256 **clP_AVX; /* handy pointers to cond likes, AVX version */ 1325 #endif 1326 #endif 1327 MrBFlt **cijks; /* space for cijks */ 1328 int **condLikeIndex; /* index to cond like space for nodes & chains */ 1329 int *condLikeScratchIndex; /* index to scratch space for node cond likes */ 1330 int **unscaledNodes; /* # unscaled nodes at a node */ 1331 int *unscaledNodesScratch; /* scratch # unscaled nodes at a node */ 1332 int **nodeScalerIndex; /* index to scaler space for nodes & chains */ 1333 int *nodeScalerScratchIndex; /* index to scratch space for node scalers */ 1334 int *siteScalerIndex; /* index to site scaler space for chains */ 1335 int siteScalerScratchIndex; /* index to scratch space for site scalers */ 1336 int **tiProbsIndex; /* index to ti probs for branches & chains */ 1337 int *tiProbsScratchIndex; /* index to scratch space for branch ti probs */ 1338 int *cijkIndex; /* index to cijks for chains */ 1339 int cijkScratchIndex; /* index to scratch space for cijks */ 1340 int numCondLikes; /* number of cond like arrays */ 1341 int numScalers; /* number of scaler arrays */ 1342 int numTiProbs; /* number of ti prob arrays */ 1343 int condLikeLength; /* length of cond like array (incl. ti cats) */ 1344 int tiProbLength; /* length of ti prob array */ 1345 MrBFlt lnLike[MAX_CHAINS]; /* log like for chain */ 1346 CLFlt *ancStateCondLikes; /* ancestral state cond like array */ 1347 1348 /* Likelihood function pointers */ 1349 LikeDownFxn CondLikeDown; /* function for calculating partials */ 1350 LikeRootFxn CondLikeRoot; /* function for calculating partials at root */ 1351 LikeScalerFxn CondLikeScaler; /* function for scaling partials */ 1352 LikeFxn Likelihood; /* function for getting cond likes for tree */ 1353 TiProbFxn TiProbs; /* function for calculating transition probs */ 1354 LikeUpFxn CondLikeUp; /* final-pass calculation of cond likes */ 1355 PrintAncStFxn PrintAncStates; /* function for sampling ancestral states */ 1356 StateCodeFxn StateCode; /* function for getting states from codes */ 1357 PrintSiteRateFxn PrintSiteRates; /* function for samling site rates */ 1358 PosSelProbsFxn PosSelProbs; /* function for sampling pos. selection probs */ 1359 SiteOmegasFxn SiteOmegas; /* function for sampling site omega values */ 1360 1361 /* Report variables */ 1362 int printAncStates; /* should ancestral states be printed (YES/NO) */ 1363 int printSiteRates; /* should site rates be printed (YES/NO) */ 1364 int printPosSel; /* should selection be printed (YES/NO) */ 1365 int printSiteOmegas; /* should site omegas be printed (YES/NO) */ 1366 1367 /* likelihood calculator flags and variables */ 1368 int useBeagle; /* use Beagle for this partition? */ 1369 int useBeagleMultiPartitions; /* use one Beagle instance for all partitions? */ 1370 int useVec; /* use SSE for this partition? */ 1371 int* rescaleFreq; /* rescale frequency for each chain */ 1372 1373 #if defined (BEAGLE_ENABLED) 1374 /* Beagle variables */ 1375 int useBeagleResource; /* try to use this BEAGLE resource number */ 1376 MrBFlt* branchLengths; /* array of branch lengths for Beagle */ 1377 MrBFlt* inRates; /* array of category rates for Beagle */ 1378 int* tiProbIndices; /* array of trans prob indices for Beagle */ 1379 MrBFlt* logLikelihoods; /* array of log likelihoods from Beagle */ 1380 int beagleInstance; /* beagle instance for division */ 1381 MrBFlt* inWeights; /* array of weights for Beagle root likelihood */ 1382 int* bufferIndices; /* array of partial indices for root likelihood */ 1383 int* eigenIndices; /* array of eigen indices for root likelihood */ 1384 int* childBufferIndices; /* array of child partial indices (unrooted) */ 1385 int* childTiProbIndices; /* array of child ti prob indices (unrooted) */ 1386 int* cumulativeScaleIndices; /* array of cumulative scale indices */ 1387 int rescaleBeagleAll; /* set to rescale all nodes */ 1388 int rescaleFreqOld; /* holds rescale frequency of current state */ 1389 int rescaleFreqNew; /* holds temporary new rescale frequency */ 1390 int recalculateScalers; /* shoud we recalculate scalers for current state YES/NO */ 1391 int* successCount; /* count of successful computations since last reset of scalers */ 1392 int** isScalerNode; /* for each node and chain set to YES if scaled node */ 1393 int* isScalerNodeScratch; /* scratch space to hold isScalerNode of proposed state*/ 1394 long* beagleComputeCount; /* count of number of calls to likelihood */ 1395 int divisionIndex; /* division index number */ 1396 BeagleOperation* operations; /* array of operations to be sent to Beagle */ 1397 int opCount; /* partial likelihood operations count */ 1398 int* scaleFactorsOps; /* array of scaler indices for Beagle operations*/ 1399 #if defined (BEAGLE_V3_ENABLED) 1400 int numCharsAll; /* number of compressed chars for all divisions */ 1401 MrBFlt* logLikelihoodsAll; /* array of log likelihoods for all divisions */ 1402 int* cijkIndicesAll; /* cijk array for all divisions */ 1403 int* categoryRateIndicesAll; /* category rate array for all divisions */ 1404 BeagleOperationByPartition* operationsAll; /* array of all operations across divisions */ 1405 BeagleOperationByPartition* operationsByPartition; /* array of division operations to be sent to Beagle */ 1406 #endif /* BEAGLE_V3_ENABLED */ 1407 #endif /* BEAGLE_ENABLED */ 1408 1409 } ModelInfo; 1410 1411 typedef struct sumt 1412 { 1413 int *absentTaxa; /* information on absent taxa */ 1414 int brlensDef; /* branch lengths defined? */ 1415 char sumtFileName[100]; /* name of input file */ 1416 char sumtOutfile[120]; /* name of output file */ 1417 char curFileName[120]; /* name of file being processed */ 1418 int burnin; /* actual burnin when parsing tree files */ 1419 char sumtConType[100]; /* consensus tree type */ 1420 int calcTreeprobs; /* should the individual tree probs be calculated*/ 1421 int showSumtTrees; /* should the individual tree probs be shown */ 1422 int numRuns; /* number of independent analyses to summarize */ 1423 int numTrees; /* number of tree params to summarize */ 1424 int orderTaxa; /* order taxa in trees? */ 1425 MrBFlt minPartFreq; /* minimum part. freq. for overall diagnostics */ 1426 int table; /* show table of partition frequencies? */ 1427 int summary; /* show summary diagnostics ? */ 1428 int showConsensus; /* show consensus trees ? */ 1429 int consensusFormat; /* format of consensus tree */ 1430 PolyTree *tree; /* for storing tree read from file */ 1431 int *order; /* for storing topology read from file */ 1432 int orderLen; /* length of order array */ 1433 int numTreesInLastBlock; /* number of trees in last block */ 1434 int numTreesEncountered; /* number of trees encounted in total */ 1435 int numTreesSampled; /* number of sampled trees in total */ 1436 int isRooted; /* is sumt tree rooted ? */ 1437 int isRelaxed; /* is sumt tree a relaxed clock tree ? */ 1438 int isClock; /* is sumt tree a clock tree ? */ 1439 int isCalibrated; /* is sumt tree calibrated ? */ 1440 int nESets; /* number of event sets */ 1441 int nBSets; /* number of branch rate sets */ 1442 char **bSetName; /* name of effective branch length sets */ 1443 char **eSetName; /* name of event sets */ 1444 int popSizeSet; /* do sumt trees have population size set? */ 1445 char *popSizeSetName; /* name of population size set */ 1446 int BitsLongsNeeded; /* number of safe longs needed for taxon bits */ 1447 int runId; /* id of run being processed */ 1448 int numTaxa; /* number of sumt taxa */ 1449 char **taxaNames; /* names of sumt taxa */ 1450 int *numFileTrees; /* number of trees per file */ 1451 int *numFileTreesSampled; /* number of trees sampled per file */ 1452 int HPD; /* use highest posterior density? */ 1453 int saveBrParams; /* save node/branch parameters? */ 1454 MrBFlt minBrParamFreq; /* threshold for printing branch params to file */ 1455 } Sumt; 1456 1457 /* formats for consensus trees */ 1458 #define SIMPLE 0 1459 #define FIGTREE 1 1460 1461 typedef struct comptree 1462 { 1463 char comptFileName1[120]; /* name of first input file */ 1464 char comptFileName2[120]; /* name of second input file */ 1465 char comptOutfile[120]; /* name of output file */ 1466 int burnin; /* actual burnin used when parsing tree files */ 1467 MrBFlt minPartFreq; /* use partitions with frequency >= minPartFreq */ 1468 } Comptree; 1469 1470 typedef struct sump 1471 { 1472 char sumpFileName[100]; /* name of input file */ 1473 char sumpOutfile[120]; /* name of output file */ 1474 //int plot; /* output plot (y/n)? */ 1475 int table; /* output table (y/n)? */ 1476 int margLike; /* output marginal likelihood (y/n)? */ 1477 int numRuns; /* number of independent analyses to summarize */ 1478 int allRuns; /* should data for all runs be printed (yes/no)? */ 1479 int HPD; /* use highest posterior density? */ 1480 MrBFlt minProb; /* cut-off for model probabilities to show */ 1481 } Sump; 1482 1483 typedef struct sumss 1484 { 1485 //int plot; /* output plot (y/n)? */ 1486 int numRuns; /* number of independent analyses to summarize */ 1487 int allRuns; /* should data for all runs be printed (yes/no)? */ 1488 int stepToPlot; /* Which step to plot in the step plot */ 1489 int askForMorePlots; /* Should user be asked to plot for different discardfraction (y/n)? */ 1490 int smoothing; /* An integer indicating number of neighbors to average over when dooing smoothing of curvs on plots */ 1491 MrBFlt discardFraction; /* Proportion of samples discarded when ploting step plot.*/ 1492 } Sumss; 1493 1494 typedef struct plot 1495 { 1496 char plotFileName[120]; /* name of input file */ 1497 char parameter[100]; /* parameter(s) to be plotted */ 1498 char match[100]; /* whether the match needs to be perfect */ 1499 } Plot; 1500 1501 typedef struct 1502 { 1503 int numTrees; /* number of trees to reassemble */ 1504 int numRuns; /* number of runs to reassemble */ 1505 } ReassembleInfo; 1506 1507 typedef struct doublet 1508 { 1509 BitsLong first, second; 1510 } Doublet; 1511 1512 typedef struct matrix 1513 { 1514 BitsLong *origin; 1515 int rowSize; 1516 int nRows; 1517 int column; 1518 int row; 1519 } Matrix; 1520 1521 typedef struct charinfo 1522 { 1523 int dType; 1524 int cType; 1525 int nStates; 1526 int constant[10]; 1527 int singleton[10]; 1528 int variable; 1529 int informative; 1530 } CharInfo; 1531 1532 typedef struct 1533 { 1534 int isExcluded; /* is the character excluded */ 1535 int numStates; /* number of observed states for the character */ 1536 int charType; /* type of character */ 1537 int isMissAmbig; /* is the character missing or ambiguous */ 1538 int ctype; /* ordering of character */ 1539 int charId; /* char ID index for doublet and codon models */ 1540 int pairsId; /* char ID for doublets */ 1541 int bigBreakAfter; /* is there a large break after this character */ 1542 } 1543 CharInformation; 1544 1545 typedef struct 1546 { 1547 int isDeleted; /* is the taxon deleted */ 1548 int charCount; /* count holder */ 1549 } 1550 TaxaInformation; 1551 1552 typedef struct 1553 { 1554 MrBFlt curScore; 1555 MrBFlt minScore; 1556 MrBFlt totalScore; 1557 MrBFlt stopScore; 1558 MrBFlt warp; 1559 TreeNode **leaf; 1560 TreeNode **vertex; 1561 } 1562 TreeInfo; 1563 1564 typedef struct 1565 { 1566 int allavailable; 1567 } 1568 ShowmovesParams; 1569 1570 typedef struct 1571 { 1572 int numNames; 1573 char **names; 1574 } 1575 NameSet; 1576 1577 /* global variables */ 1578 extern int abortMove; /* flag determining whether to abort move */ 1579 extern int *activeParams[NUM_LINKED]; /* a table holding the parameter status */ 1580 extern int *activeParts; /* partitions changes should apply to */ 1581 extern int autoClose; /* autoclose */ 1582 extern int autoOverwrite; /* Overwrite or append outputfiles when nowarnings=yes */ 1583 extern int chainHasAdgamma; /* indicates if chain has adgamma HMMs */ 1584 extern Chain chainParams; /* holds parameters for Markov chain */ 1585 extern CharInformation *charInfo; /* holds critical information about characters */ 1586 extern char **charSetNames; /* holds names of character sets */ 1587 extern int *compCharPos; /* char position in compressed matrix */ 1588 extern int *compColPos; /* column position in compressed matrix */ 1589 extern BitsLong *compMatrix; /* compressed character matrix */ 1590 extern int compMatrixRowSize; /* row size of compressed matrix */ 1591 extern Comptree comptreeParams; /* holds parameters for comparetree command */ 1592 extern char **constraintNames; /* holds names of constraints */ 1593 extern BitsLong **definedConstraint; /* holds information about defined constraints */ 1594 extern BitsLong **definedConstraintTwo; /* bitfields representing second taxa sets of defined constraints (for PARTIAL constraints) */ 1595 extern BitsLong **definedConstraintPruned; /* bitfields representing taxa sets of defined constraints after delited taxa are removed */ 1596 extern BitsLong **definedConstraintTwoPruned; /* bitfields representing second taxa sets of defined constraints after delited taxa are removed (for PARTIAL constraints) */ 1597 extern int dataType; /* type of data */ 1598 extern Calibration defaultCalibration; /* default model settings */ 1599 extern ModelParams defaultModel; /* default model settings */ 1600 extern int defChars; /* flag for whether number of characters is known*/ 1601 extern int defMatrix; /* flag for whether matrix is successfull read */ 1602 extern int defPairs; /* flag for whether constraints on tree are read */ 1603 extern int defPartition; /* flag for whether character partition is read */ 1604 extern int defTaxa; /* are taxon labels defined ? */ 1605 extern Doublet doublet[16]; /* holds information on states for doublets */ 1606 extern int echoMB; /* flag used by Manual to prevent echoing */ 1607 extern BitsLong expecting; /* variable denoting expected token type */ 1608 extern int fileNameChanged; /* has file name been changed? */ 1609 extern int foundNewLine; /* whether a new line has been found */ 1610 extern char gapId; /* gap character Id */ 1611 extern RandLong globalSeed; /* seed that is initialized at start up */ 1612 extern char **headerNames; /* string to hold headers in sump and plot */ 1613 extern int inComment; /* flag for whether input stream is commented */ 1614 extern int inferAncStates; /* should ancestral states be inferred (y/n) */ 1615 extern int inferSiteOmegas; /* should site omega values be inferred (y/n) */ 1616 extern int inferSiteRates; /* should site rates be inferred (y/n) */ 1617 extern int inferPosSel; /* should positive selection be inferred (y/n) */ 1618 extern char inputFileName[100]; /* input (NEXUS) file name */ 1619 extern int inTreesBlock; /* are we in the sumt block */ 1620 extern int inValidCommand; /* a useful flag set whenever you enter a cmd */ 1621 extern int isInAmbig, isInPoly; /* flags whether we are within () or {} */ 1622 extern int isMixed; /* flags whether dataset is mixed */ 1623 extern int inMrbayesBlock; /* flag for whether we are in a mrbayes block */ 1624 extern int *intValues; /* integer values of parameters */ 1625 extern int isTaxsetDef; /* is a taxon set defined */ 1626 extern int isTranslateDef; /* is a translation block defined */ 1627 extern int isTranslateDiff; /* is translate different from current taxaset? */ 1628 extern int *linkTable[NUM_LINKED]; /* how parameters are linked across parts */ 1629 extern int localOutGroup; /* outgroup for non-excluded taxa */ 1630 extern char **localTaxonNames; /* points to names of non-excluded taxa */ 1631 extern FILE *logFileFp; /* file pointer to log file */ 1632 extern char logFileName[100]; /* name of the log file */ 1633 extern int logToFile; /* should screen output be logged to a file */ 1634 extern char manFileName[100]; /* name of man file */ 1635 extern char matchId; /* mach character Id */ 1636 extern int *matrix; /* matrix containing original data */ 1637 extern int matrixHasPoly; /* flag for whether matrix has polymorphisms */ 1638 extern int memAllocs[NUM_ALLOCS]; /* allocated memory flags */ 1639 extern int mode; /* mode of program (interactive/noninteractive) */ 1640 extern char **modelIndicatorParams; /* model indicator params */ 1641 extern char ***modelElementNames; /* names for component models */ 1642 extern MCMCMove **moves; /* vector of applicable moves */ 1643 extern MoveType moveTypes[NUM_MOVE_TYPES]; /* holds information on the move types */ 1644 extern char missingId; /* missing character Id */ 1645 extern Tree **mcmcTree; /* pointers to mcmc trees */ 1646 extern Model *modelParams; /* holds model params for partitions */ 1647 extern ModelInfo *modelSettings; /* stores important info on model params */ 1648 extern int nBitsInALong; /* number of bits in a BitsLong */ 1649 extern Calibration *nodeCalibration; /* holds information about node calibrations */ 1650 extern int noWarn; /* no warnings on overwriting files */ 1651 extern int nPThreads; /* number of pthreads to use */ 1652 extern int numActiveLocks; /* number of active, locked nodes */ 1653 extern int numApplicableMoves; /* number of moves applicable to parameters */ 1654 extern int numChar; /* number of characters in character matrix */ 1655 extern int numCharSets; /* holds number of character sets */ 1656 extern int numComments; /* number of nested comments */ 1657 extern int numCompressedChars; /* number of compressed characters */ 1658 extern int numCurrentDivisions; /* number of partitions of data */ 1659 extern int numDefinedConstraints; /* number of constraints defined */ 1660 extern enum ConstraintType *definedConstraintsType; /* Store type of constraint */ 1661 extern int numDefinedPartitions; /* number of partitions defined */ 1662 extern int numDefinedSpeciespartitions; /* number of species partitions defined */ 1663 extern int numGlobalChains; /* number of global chains */ 1664 extern int numLocalTaxa; /* number of non-excluded taxa */ 1665 extern int numLocalChar; /* number of non-excluded characters */ 1666 extern int numMoveTypes; /* the number of move types */ 1667 extern int numOpenExeFiles; /* number of execute files open */ 1668 extern int numParams; /* number of parameters in model */ 1669 extern int numDivisions; /* number of current divisions */ 1670 extern int numPrintParams; /* number of substitution model parameters to print */ 1671 extern int numPrintTreeParams; /* number of tree model parameters to print */ 1672 extern CLFlt *numSitesOfPat; /* no. sites of each pattern */ 1673 extern int numSpecies; /* number of species in current speciespartition */ 1674 extern int numTaxa; /* number of taxa in character matrix */ 1675 extern int numTaxaSets; /* holds number of taxa sets */ 1676 extern int numTopologies; /* number of topologies for one chain and state */ 1677 extern int numTranslates; /* number of taxa in active translate block */ 1678 extern int numTrees; /* number of trees for one chain and state */ 1679 extern int numUserTrees; /* number of defined user trees */ 1680 extern int *numVars; /* number of variables in setting arrays */ 1681 extern int *origChar; /* index from compressed char to original char */ 1682 extern int outGroupNum; /* number of outgroup taxon */ 1683 extern ParmInfo paramTable[]; /* information on parameters */ 1684 extern MrBFlt *paramValues; /* values of parameters */ 1685 extern int **partitionId; /* holds information about defined partitions */ 1686 extern char **partitionNames; /* hold names of partitions (first is "default") */ 1687 extern MrBFlt *parameterValues; /* vector holding sump or plot parameters */ 1688 extern Param *params; /* vector of parameters in model */ 1689 extern int partitionNum; /* index of current partition */ 1690 extern Plot plotParams; /* holds parameters for plot command */ 1691 extern int precision; /* precision of samples and summary stats */ 1692 extern int *printAncStates; /* divisions to print anc states for */ 1693 extern int quitOnError; /* quit on error? */ 1694 extern int readComment; /* should we read comment (looking for &)? */ 1695 extern int readWord; /* should we read a word next? */ 1696 extern ReassembleInfo reassembleParams; /* holds parameters for reassemble command */ 1697 extern int replaceLogFile; /* should logfile be replace/appended to */ 1698 extern RandLong runIDSeed; /* seed used only for generating run ID [stamp] */ 1699 extern BitsLong bitsLongWithAllBitsSet; /* a BitsLong with all bits set, for bit ops */ 1700 extern int setUpAnalysisSuccess; /* Set to YES if analysis is set without error */ 1701 extern int scientific; /* use scientific format for samples ? */ 1702 extern ShowmovesParams showmovesParams; /* holds parameters for Showmoves command */ 1703 extern char spacer[10]; /* holds blanks for printing indentations */ 1704 extern NameSet *speciesNameSets; /* hold species name sets, one for each speciespartition */ 1705 extern int **speciespartitionId; /* holds info about defined speciespartitions */ 1706 extern char **speciespartitionNames; /* hold names of speciespartitions (first is "default") */ 1707 extern int speciespartitionNum; /* index of current species partition */ 1708 extern char stamp[11]; /* holds a unique identifier for each analysis */ 1709 extern RandLong swapSeed; /* seed used only for determining which to swap */ 1710 extern int state[MAX_CHAINS]; /* state of chain */ 1711 extern MrBFlt *stdStateFreqs; /* std char state frequencies */ 1712 extern int *stdType; /* compressed std char type: ord, unord, irrev */ 1713 extern Sump sumpParams; /* holds parameters for sump command */ 1714 extern char sumpToken[]; /* string holding a .p file token */ 1715 extern char *sumpTokenP; /* pointer to a .p file token */ 1716 extern Sumt sumtParams; /* holds parameters for sumt command */ 1717 extern Sumss sumssParams; /* holds parameters for sumss command */ 1718 extern int stdStateFreqsRowSize; /* row size for stdStateFreqs */ 1719 extern int *sympiIndex; /* sympi state freq index for multistate chars */ 1720 extern TaxaInformation *taxaInfo; /* holds critical information about taxa */ 1721 extern char **taxaNames; /* holds name of taxa */ 1722 extern char **taxaSetNames; /* holds names of taxa sets */ 1723 extern BitsLong **taxaSet; /* holds information about defined taxasets */ 1724 extern int *tempActiveConstraints; /* info on the active constraints in prset */ 1725 extern int *tempLinkUnlink[NUM_LINKED]; /* for changing parameter linkage */ 1726 extern int *tempLinkUnlinkVec; /* for changing parameter linkage */ 1727 extern MrBFlt *tempNum; /* vector of numbers used for setting arrays */ 1728 extern int *tempSet; /* temporarily holds defined character set */ 1729 extern int theAmbigChar; /* int containing ambiguous character */ 1730 extern int *tiIndex; /* compressed std char ti index */ 1731 extern Calibration *tipCalibration; /* holds tip calibrations */ 1732 extern char **transFrom; /* translation block information */ 1733 extern char **transTo; /* translation block information */ 1734 extern int userBrlensDef; /* are the branch lengths on user tree defined */ 1735 extern int userLevel; /* the level of the user */ 1736 extern PolyTree *userTree[]; /* array of user trees */ 1737 extern char workingDir[100]; /* working directory */ 1738 #if defined (BEAGLE_ENABLED) 1739 extern int tryToUseBEAGLE; /* try to use the BEAGLE library */ 1740 extern long beagleFlags; /* BEAGLE requirement flags */ 1741 extern int beagleResourceNumber; /* BEAGLE resource number */ 1742 extern int* beagleResource; /* BEAGLE resource list */ 1743 extern int beagleResourceCount; /* BEAGLE resource list length */ 1744 extern int beagleInstanceCount; /* total number of BEAGLE instances */ 1745 extern int beagleScalingScheme; /* BEAGLE dynamic scaling */ 1746 extern int beagleScalingFrequency; /* BEAGLE rescaling frequency */ 1747 extern int recalcScalers; /* shoud we recalculate scalers for one of divisions for current state YES/NO */ 1748 #if defined (BEAGLE_V3_ENABLED) 1749 extern int beagleThreadCount; /* max number of BEAGLE CPU threads */ 1750 extern int beagleAllFloatTips; /* use floating-point represantion for all tips */ 1751 #endif 1752 #endif 1753 1754 /* Aamodel parameters */ 1755 extern MrBFlt aaJones[20][20]; /* rates for Jones model */ 1756 extern MrBFlt aaDayhoff[20][20]; /* rates for Dayhoff model */ 1757 extern MrBFlt aaMtrev24[20][20]; /* rates for mtrev24 model */ 1758 extern MrBFlt aaMtmam[20][20]; /* rates for mtmam model */ 1759 extern MrBFlt aartREV[20][20]; /* rates for rtREV model */ 1760 extern MrBFlt aaWAG[20][20]; /* rates for WAG model */ 1761 extern MrBFlt aacpREV[20][20]; /* rates for aacpREV model */ 1762 extern MrBFlt aaVt[20][20]; /* rates for VT model */ 1763 extern MrBFlt aaBlosum[20][20]; /* rates for Blosum62 model */ 1764 extern MrBFlt aaLG[20][20]; /* rates for LG model */ 1765 extern MrBFlt jonesPi[20]; /* stationary frequencies for Jones model */ 1766 extern MrBFlt dayhoffPi[20]; /* stationary frequencies for Dayhoff model */ 1767 extern MrBFlt mtrev24Pi[20]; /* stationary frequencies for mtrev24 model */ 1768 extern MrBFlt mtmamPi[20]; /* stationary frequencies for mtmam model */ 1769 extern MrBFlt rtrevPi[20]; /* stationary frequencies for rtREV model */ 1770 extern MrBFlt wagPi[20]; /* stationary frequencies for WAG model */ 1771 extern MrBFlt cprevPi[20]; /* stationary frequencies for aacpREV model */ 1772 extern MrBFlt vtPi[20]; /* stationary frequencies for VT model */ 1773 extern MrBFlt blosPi[20]; /* stationary frequencies for Blosum62 model */ 1774 extern MrBFlt lgPi[20]; /* stationary frequencies for LG model */ 1775 1776 #if defined (PRINT_DUMP) 1777 FILE *dumpFile; /* for debugging logs */ 1778 #endif 1779 1780 #if defined (MPI_ENABLED) 1781 extern int proc_id; /* process ID (0, 1, ..., num_procs-1) */ 1782 extern int num_procs; /* number of active processors */ 1783 extern MrBFlt myStateInfo[7]; /* likelihood/prior/heat/ran/moveInfo vals of me */ 1784 extern MrBFlt partnerStateInfo[7]; /* likelihood/prior/heat/ran/moveInfo vals of partner */ 1785 #endif 1786 1787 #endif /* __BAYES_H__ */ 1788