src/bayes.h

   1 #ifndef __BAYES_H__
   2 #define __BAYES_H__
   3
   4 #include <assert.h>
   5 #include <ctype.h>
   6 #include <float.h>
   7 #include <limits.h>
   8 #include <math.h>
   9 #include <memory.h>
  10 #include <stdio.h>
  11 #include <stdlib.h>
  12 #include <string.h>
  13 #include <stdarg.h>
  14 #include <time.h>
  15
  16 #ifdef USECONFIG_H
  17 #   include "config.h"
  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
  24 #   undef  MPI_ENABLED
  25 #   undef  BEAGLE_ENABLED
  26 #   undef  FAST_LOG
  27 #endif
  28
  29 #if defined (MPI_ENABLED)
  30 #include "mpi.h"
  31 #endif
  32
  33 #if defined (BEAGLE_ENABLED)
  34 #include "libhmsbeagle/beagle.h"
  35 #endif
  36
  37 /* uncomment the following line when releasing, also modify the VERSION_NUMBER below */
  38 #define RELEASE
  39 #ifdef RELEASE
  40 #define VERSION_NUMBER  "3.2.6"
  41 #else
  42 #define VERSION_NUMBER  "3.2.7-svn"
  43 #endif
  44
  45 #if !defined (UNIX_VERSION) && !defined (WIN_VERSION) && !defined (MAC_VERSION)
  46 #  ifdef __MWERKS__
  47 #    define MAC_VERSION
  48 #  elif defined __APPLE__
  49 #    define MAC_VERSION
  50 #  else
  51 #    define WIN_VERSION
  52 #  endif
  53 #endif
  54
  55 #if defined (WIN_VERSION)
  56 #   include <windows.h>
  57 #   include <winbase.h>
  58 #endif
  59
  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.
  65
  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
  70    some systems.
  71
  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.
  76
  77    FR 2013-07-06
  78  */
  79 typedef unsigned long BitsLong;
  80 typedef long RandLong;
  81
  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 */
  90
  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)
  95 #    define MS_VCPP_SSE
  96 #    include <xmmintrin.h>
  97 #  else
  98 #    define GCC_SSE
  99 #    undef ICC_SSE
 100 #    include <xmmintrin.h>
 101 #  endif
 102 #endif
 103
 104 /* For comparing floating points: two values are the same if the absolute difference is less then
 105    this value.
 106 */
 107 #ifndef ETA
 108 #define ETA (1E-30)
 109 #endif
 110
 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);
 114 #endif
 115 */
 116
 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.
 121 */
 122 #ifdef __sgi
 123 #define TEMPSTRSIZE 1000
 124 #else
 125 #define TEMPSTRSIZE 200
 126 #endif
 127
 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 */
 130 #undef  NO_ERROR
 131 #undef  ERROR
 132 #define NO_ERROR                0
 133 #define ERROR                   1
 134 #define NO_ERROR_QUIT           2
 135 #define ABORT                   3
 136 #define SKIP_COMMAND            4
 137
 138 #undef  FALSE
 139 #undef  TRUE
 140 #define FALSE                   0
 141 #define TRUE                    1
 142
 143 #define NO                      0
 144 #define YES                     1
 145
 146 #define UP                      0
 147 #define DOWN                    1
 148
 149 #define UPPER                   0
 150 #define MIDDLE                  1
 151 #define LOWER                   2
 152
 153 #define NONINTERACTIVE          0
 154 #define INTERACTIVE             1
 155
 156 #define STANDARD_USER           1
 157 #define DEVELOPER               3
 158
 159 #define DIFFERENT               0
 160 #define SAME                    1
 161 #define CONSISTENT_WITH         2
 162
 163 #define LINETERM_UNIX           0
 164 #define LINETERM_MAC            1
 165 #define LINETERM_DOS            2
 166
 167 #define SCREENWIDTH             60
 168 #define SCREENWIDTH2            61
 169
 170 #define AVGSTDDEV               0
 171 #define MAXSTDDEV               1
 172
 173 #define NONE                    0
 174 #define DNA                     1
 175 #define RNA                     2
 176 #define PROTEIN                 3
 177 #define RESTRICTION             4
 178 #define STANDARD                5
 179 #define MIXED                   6
 180 #define CONTINUOUS              7
 181
 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
 190 #define AAMODEL_VT              8
 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 */
 195
 196 #define NUCMODEL_4BY4           0
 197 #define NUCMODEL_DOUBLET        1
 198 #define NUCMODEL_CODON          2
 199 #define NUCMODEL_AA             3
 200
 201 #define NST_MIXED              -1  /* anything other than 1, 2, or 6 */
 202
 203 #define MISSING                 10000000
 204 #define GAP                     10000001
 205
 206 #define UNORD                   0
 207 #define ORD                     1
 208 #define DOLLO                   2
 209 #define IRREV                   3
 210
 211 #define IN_CMD                  0
 212 #define IN_FILE                 1
 213
 214 #define NOTHING                 0
 215 #define COMMAND                 1
 216 #define PARAMETER               2
 217 #define EQUALSIGN               3
 218 #define COLON                   4
 219 #define SEMICOLON               5
 220 #define COMMA                   6
 221 #define POUNDSIGN               7
 222 #define QUESTIONMARK            8
 223 #define DASH                    9
 224 #define LEFTPAR                 10
 225 #define RIGHTPAR                11
 226 #define LEFTCOMMENT             12
 227 #define RIGHTCOMMENT            13
 228 #define ALPHA                   14
 229 #define NUMBER                  15
 230 #define RETURNSYMBOL            16
 231 #define ASTERISK                17
 232 #define BACKSLASH               18
 233 #define FORWARDSLASH            19
 234 #define EXCLAMATIONMARK         20
 235 #define PERCENT                 21
 236 #define QUOTATIONMARK           22
 237 #define WEIRD                   23
 238 #define UNKNOWN_TOKEN_TYPE      24
 239 #define LEFTCURL                25
 240 #define RIGHTCURL               26
 241 #define DOLLAR                  27
 242 #define AMPERSAND               28
 243 #define VERTICALBAR             29
 244
 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
 288
 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
 295
 296 #define CMD_STRING_LENGTH       100000
 297
 298 #define pos(i,j,n)              ((i)*(n)+(j))
 299
 300 #define NUM_ALLOCS               100
 301
 302 #define ALLOC_MATRIX             0
 303 #define ALLOC_CHARINFO           2
 304 #define ALLOC_CHARSETS           3
 305 #define ALLOC_TAXA               4
 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
 343 #define ALLOC_PBF                68
 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
 358 #define ALLOC_SS                 90
 359 #define ALLOC_SAMPLEFOSSILSLICE  91
 360
 361 #define LINKED                  0
 362 #define UNLINKED                1
 363
 364 /*paramType*/
 365 #define NUM_LINKED              31
 366 #define P_TRATIO                0
 367 #define P_REVMAT                1
 368 #define P_OMEGA                 2
 369 #define P_PI                    3
 370 #define P_SHAPE                 4
 371 #define P_PINVAR                5
 372 #define P_CORREL                6
 373 #define P_SWITCH                7
 374 #define P_RATEMULT              8
 375 #define P_TOPOLOGY              9
 376 #define P_BRLENS                10
 377 #define P_SPECRATE              11
 378 #define P_EXTRATE               12
 379 #define P_FOSLRATE              13
 380 #define P_POPSIZE               14
 381 #define P_AAMODEL               15
 382 #define P_BRCORR                16
 383 #define P_BRSIGMA               17
 384 #define P_GROWTH                18
 385 #define P_CPPMULTDEV            19
 386 #define P_CPPRATE               20
 387 #define P_CPPEVENTS             21
 388 #define P_TK02VAR               22
 389 #define P_TK02BRANCHRATES       23
 390 #define P_IGRVAR                24
 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 */
 398
 399 // #define CPPm                 0       /* CPP rate multipliers */
 400 // #define CPPi                 1       /* CPP independent rates */
 401 #define RCL_TK02                0
 402 #define RCL_IGR                 1       /* type of mixed relaxed clock model */
 403
 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] */
 406
 407 // #define PARAM_NAME_SIZE      400
 408
 409 typedef void * VoidPtr;
 410 typedef int (*CmdFxn)(void);
 411 typedef int (*ParmFxn)(char *, char *);
 412
 413 /* typedef for a ln prior prob fxn */
 414 typedef MrBFlt (*LnPriorProbFxn)(MrBFlt val, MrBFlt *priorParams);
 415
 416 /* typedef for a ln prior prob ratio fxn */
 417 typedef MrBFlt (*LnPriorRatioFxn)(MrBFlt newVal, MrBFlt oldVal, MrBFlt *priorParams);
 418
 419 typedef struct
 420     {
 421     MrBFlt          sum;            /* sum of standard deviations */
 422     MrBFlt          max;            /* maximum standard deviation */
 423     MrBFlt          numPartitions;
 424     MrBFlt          numSamples;
 425     MrBFlt          avgStdDev;
 426     MrBFlt          **pair;
 427     } STATS;
 428
 429 /* enumeration for calibration prior type */
 430 enum CALPRIOR
 431     {
 432     unconstrained,
 433     fixed,
 434     uniform,
 435     offsetExponential,
 436     truncatedNormal,
 437     logNormal,
 438     offsetLogNormal,
 439     standardGamma,
 440     offsetGamma
 441     };
 442
 443 enum ConstraintType
 444     {
 445     PARTIAL,
 446     NEGATIVE,
 447     HARD
 448     };
 449
 450 enum CodingType
 451     {
 452     ALL                 = 0,
 453     NOABSENCESITES      = 1,
 454     NOPRESENCESITES     = 2,
 455     VARIABLE            = 3,
 456     NOSINGLETONPRESENCE = 4,
 457     NOSINGLETONABSENCE  = 8,
 458     NOSINGLETONS        = 12,
 459     INFORMATIVE         = 15
 460     };
 461
 462 /* typedef for calibration */
 463 typedef struct calibration
 464     {
 465     char                name[100];
 466     enum CALPRIOR       prior;
 467     MrBFlt              priorParams[3];
 468     LnPriorProbFxn      LnPriorProb;
 469     LnPriorRatioFxn     LnPriorRatio;
 470     MrBFlt              min;
 471     MrBFlt              max;
 472     }
 473     Calibration;
 474
 475 /* typedef for tree (topology) list element */
 476 typedef struct element
 477     {
 478     struct element *next;
 479     int             *order;
 480     } TreeListElement;
 481
 482 /* typedef for list of trees (topologies) */
 483 typedef struct
 484     {
 485     TreeListElement *first;
 486     TreeListElement *last;
 487     } TreeList;
 488
 489 /* typedef for packed tree */
 490 typedef struct
 491     {
 492     int     *order;
 493     MrBFlt  *brlens;
 494     } PackedTree;
 495
 496 /* typedef for binary tree node */
 497 /* NOTE: Any variable added here must also be copied in CopyTrees */
 498 typedef struct node
 499     {
 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                */
 517     }
 518     TreeNode;
 519
 520 /* typedef for binary tree */
 521 typedef struct
 522     {
 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 */
 543     }
 544     Tree;
 545
 546 /* typedef for node in polytomous tree */
 547 typedef struct pNode
 548     {
 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                    */
 564     }
 565     PolyNode;
 566
 567 /* typedef for polytomous tree */
 568 typedef struct
 569     {
 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             */
 596     }
 597     PolyTree;
 598
 599 /* struct for holding model parameter info for the mcmc run */
 600 typedef struct param
 601     {
 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                   */
 642     } Param;
 643
 644 #if defined(THREADS_ENABLED)
 645 #include <pthread.h>
 646
 647 typedef struct s_launch_struct
 648     {
 649     int chain;
 650     int division;
 651     MrBFlt* lnL;
 652     } LaunchStruct;
 653 #endif
 654
 655 /* parameter ID values */
 656 /* identifies unique model parameter x prior combinations */
 657 #define TRATIO_DIR                      1
 658 #define TRATIO_FIX                      2
 659 #define REVMAT_DIR                      3
 660 #define REVMAT_FIX                      4
 661 #define OMEGA_DIR                       5
 662 #define OMEGA_FIX                       6
 663 #define SYMPI_UNI                       7
 664 #define SYMPI_UNI_MS                    8
 665 #define SYMPI_EXP                       9
 666 #define SYMPI_EXP_MS                    10
 667 #define SYMPI_FIX                       11
 668 #define SYMPI_FIX_MS                    12
 669 #define SYMPI_EQUAL                     13
 670 #define PI_DIR                          14
 671 #define PI_USER                         15
 672 #define PI_EMPIRICAL                    16
 673 #define PI_EQUAL                        17
 674 #define PI_FIXED                        18
 675 #define SHAPE_UNI                       19
 676 #define SHAPE_EXP                       20
 677 #define SHAPE_FIX                       21
 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
 736 #define OMEGA_BUD                       80
 737 #define OMEGA_BUF                       81
 738 #define OMEGA_BED                       82
 739 #define OMEGA_BEF                       83
 740 #define OMEGA_BFD                       84
 741 #define OMEGA_BFF                       85
 742 #define OMEGA_FUD                       86
 743 #define OMEGA_FUF                       87
 744 #define OMEGA_FED                       88
 745 #define OMEGA_FEF                       89
 746 #define OMEGA_FFD                       90
 747 #define OMEGA_FFF                       91
 748 #define OMEGA_ED                        92
 749 #define OMEGA_EF                        93
 750 #define OMEGA_FD                        94
 751 #define OMEGA_FF                        95
 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
 803
 804 #if defined (BEAGLE_ENABLED)
 805 #define MB_BEAGLE_SCALE_ALWAYS          0
 806 #define MB_BEAGLE_SCALE_DYNAMIC         1
 807 #if defined (_DEBUG)
 808 #define MB_PRINT_DYNAMIC_RESCALE_FAIL_STAT
 809 #endif
 810 #endif
 811
 812 /* typedef for a MoveFxn */
 813 typedef int (MoveFxn)(Param *param, int chain, RandLong *seed, MrBFlt *lnPriorRatio, MrBFlt *lnProposalRatio, MrBFlt *mvp);
 814
 815 /* typedef for an ApplicFxn */
 816 typedef int (ApplicFxn)(Param *param);
 817
 818 /* typedef for an AutotuneFxn */
 819 typedef void (AutotuneFxn)(MrBFlt acceptanceRate, MrBFlt targetRate, int batch, MrBFlt *tuningParameter, MrBFlt minTuning, MrBFlt maxTuning);
 820
 821 /* struct holding info on each move type that the program handles */
 822 typedef struct
 823     {
 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 */
 843     } MoveType;
 844
 845 /* max number of move types */
 846 #define NUM_MOVE_TYPES 100
 847
 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 */
 852 typedef struct
 853     {
 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               */
 868     } MCMCMove;
 869
 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);
 881
 882 typedef struct cmdtyp
 883     {
 884     int         cmdNumber;
 885     char        *string;
 886     int         specialCmd;
 887     CmdFxn      cmdFxnPtr;
 888     short       numParms;
 889     short       parmList[50];
 890     int         expect;
 891     char        *cmdDescription;
 892     int         cmdUse;
 893     int         hiding;
 894     }
 895     CmdType;
 896
 897 typedef struct parm
 898     {
 899     char        *string;    /* parameter name */
 900     char        *valueList; /* list of values that could be input */
 901     ParmFxn     fp;         /* function pointer */
 902     }
 903     ParmInfo, *ParmInfoPtr;
 904
 905 typedef struct model
 906     {
 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        */
 912
 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 */
 925
 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                     */
 931
 932     char        tRatioPr[100];     /* prior for ti/tv rate ratio                   */
 933     MrBFlt      tRatioFix;
 934     MrBFlt      tRatioDir[2];
 935     char        revMatPr[100];     /* prior for GTR model                          */
 936     MrBFlt      revMatFix[6];
 937     MrBFlt      revMatDir[6];
 938     MrBFlt      revMatSymDir;      /* prior for mixed GTR subspace model           */
 939     char        aaModelPr[100];    /* prior for amino acid model                   */
 940     char        aaModel[100];
 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                              */
 946     MrBFlt      omegaFix;
 947     MrBFlt      omegaDir[2];
 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];
 959     MrBFlt      m10betaExp;
 960     MrBFlt      m10betaUni[2];
 961     MrBFlt      m10betaFix[2];
 962     MrBFlt      m10gammaExp;
 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];
 972     int         numDirParams;
 973     char        shapePr[100];      /* prior for gamma/lnorm shape parameter        */
 974     MrBFlt      shapeFix;
 975     MrBFlt      shapeUni[2];
 976     MrBFlt      shapeExp;
 977     char        pInvarPr[100];     /* prior for proportion of invariable sites     */
 978     MrBFlt      pInvarFix;
 979     MrBFlt      pInvarUni[2];
 980     char        adGammaCorPr[100]; /* prior for correlation param of adGamma model */
 981     MrBFlt      corrFix;
 982     MrBFlt      corrUni[2];
 983     char        covSwitchPr[100];  /* prior for switching rates of covarion model  */
 984     MrBFlt      covswitchFix[2];
 985     MrBFlt      covswitchUni[2];
 986     MrBFlt      covswitchExp;
 987     char        symPiPr[100];      /* prior for pi when unidentifiable states used */
 988     MrBFlt      symBetaFix;
 989     MrBFlt      symBetaUni[2];
 990     MrBFlt      symBetaExp;
 991     char        ratePr[100];       /* prior on rate for a partition                */
 992     MrBFlt      ratePrDir;
 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      */
 996     MrBFlt      brownCorrFix;
 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;
1003
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;
1008     int         numActiveLocks;
1009     char        brlensPr[100];     /* prior on branch lengths                      */
1010     int         brlensFix;         /* user tree index for fixed brlens             */
1011     MrBFlt      brlensUni[2];
1012     MrBFlt      brlensExp;
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                      */
1044     MrBFlt      popSizeFix;
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              */
1051     MrBFlt      growthFix;
1052     MrBFlt      growthUni[2];
1053     MrBFlt      growthExp;
1054     MrBFlt      growthNorm[2];
1055     char        cppRatePr[100];        /* prior on CPP rate                             */
1056     MrBFlt      cppRateFix;
1057     MrBFlt      cppRateExp;
1058     char        cppMultDevPr[100];     /* prior on CPP rate multiplier Lognormal variance */
1059     MrBFlt      cppMultDevFix;
1060     char        tk02varPr[100];        /* prior on TK02 lognormal rate variance         */
1061     MrBFlt      tk02varFix;
1062     MrBFlt      tk02varUni[2];
1063     MrBFlt      tk02varExp;
1064     char        igrvarPr[100];         /* prior on IGR gamma distribution variance      */
1065     MrBFlt      igrvarFix;
1066     MrBFlt      igrvarUni[2];
1067     MrBFlt      igrvarExp;
1068     char        mixedvarPr[100];       /* prior on mixed relaxed clock rate variance    */
1069     MrBFlt      mixedvarFix;
1070     MrBFlt      mixedvarUni[2];
1071     MrBFlt      mixedvarExp;
1072
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;
1082
1083 typedef struct chain
1084     {
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                          */
1140     } Chain;
1141
1142 typedef struct modelinfo
1143     {
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           */
1150
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     */
1158
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              */
1189
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)        */
1203
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      */
1208
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  */
1216
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 */
1222
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             */
1229
1230     /* Miscellaneous parameters */
1231     int         mark;                       /* scratch parameter                            */
1232     int         parsimonyBasedMove;         /* is parsimony-based move used (YES/NO)        */
1233
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   */
1249 #endif
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              */
1270
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      */
1283
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)       */
1289
1290     /* likelihood calculator flags */
1291     int         useBeagle;                  /* use Beagle for this partition?               */
1292     int         useSSE;                     /* use SSE for this partition?                  */
1293
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       */
1316 #endif
1317
1318     } ModelInfo;
1319
1320 typedef struct sumt
1321     {
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?                */
1364     } Sumt;
1365
1366 /* formats for consensus trees */
1367 #define SIMPLE      0
1368 #define FIGTREE     1
1369
1370 typedef struct comptree
1371     {
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  */
1377     } Comptree;
1378
1379 typedef struct sump
1380     {
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       */
1390     } Sump;
1391
1392 typedef struct sumss
1393     {
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.*/
1401     } Sumss;
1402
1403 typedef struct plot
1404     {
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         */
1408     } Plot;
1409
1410 typedef struct
1411     {
1412     int         numTrees;              /* number of trees to reassemble                 */
1413     int         numRuns;               /* number of runs to reassemble                  */
1414     } ReassembleInfo;
1415
1416 typedef struct doublet
1417     {
1418     BitsLong    first, second;
1419     } Doublet;
1420
1421 typedef struct matrix
1422     {
1423     BitsLong *origin;
1424     int rowSize;
1425     int nRows;
1426     int column;
1427     int row;
1428     } Matrix;
1429
1430 typedef struct charinfo
1431     {
1432     int dType;
1433     int cType;
1434     int nStates;
1435     int constant[10];
1436     int singleton[10];
1437     int variable;
1438     int informative;
1439     } CharInfo;
1440
1441 typedef struct
1442     {
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   */
1451     }
1452     CharInformation;
1453
1454 typedef struct
1455     {
1456     int         isDeleted;             /* is the taxon deleted                          */
1457     int         charCount;             /* count holder                                  */
1458     }
1459     TaxaInformation;
1460
1461 typedef struct
1462     {
1463     MrBFlt      curScore;
1464     MrBFlt      minScore;
1465     MrBFlt      totalScore;
1466     MrBFlt      stopScore;
1467     MrBFlt      warp;
1468     TreeNode    **leaf;
1469     TreeNode    **vertex;
1470     }
1471     TreeInfo;
1472
1473 typedef struct
1474     {
1475     int     allavailable;
1476     }
1477     ShowmovesParams;
1478
1479 typedef struct
1480     {
1481     int     numNames;
1482     char    **names;
1483     }
1484     NameSet;
1485
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 */
1657 #endif
1658 #if defined (THREADS_ENABLED)
1659 extern int              tryToUseThreads;                        /* try to use pthreads with BEAGLE library       */
1660 #endif
1661
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          */
1683
1684 #if defined (PRINT_DUMP)
1685 FILE                    *dumpFile;                   /* for debugging logs */
1686 #endif
1687
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         */
1693 #endif
1694
1695 #if defined (FAST_LOG)
1696 extern CLFlt            scalerValue[];
1697 extern CLFlt            logValue[];
1698 #endif
1699
1700 #endif  /* __BAYES_H__ */