| 1 | #ifndef __MB_H__ |
|---|
| 2 | #define __MB_H__ |
|---|
| 3 | |
|---|
| 4 | #include <stdio.h> |
|---|
| 5 | #include <float.h> |
|---|
| 6 | |
|---|
| 7 | #ifdef USECONFIG_H |
|---|
| 8 | # include "config.h" |
|---|
| 9 | #else /* some defaults that would otherwise be guessed by configure */ |
|---|
| 10 | # define PACKAGE_NAME "mrbayes" |
|---|
| 11 | # define PACKAGE_VERSION "3.2" |
|---|
| 12 | # undef HAVE_LIBREADLINE |
|---|
| 13 | # define UNIX_VERSION 1 |
|---|
| 14 | # undef FAST_LOG |
|---|
| 15 | #if !defined(XCODE_VERSION) |
|---|
| 16 | # undef _64BIT |
|---|
| 17 | #endif |
|---|
| 18 | #endif |
|---|
| 19 | |
|---|
| 20 | #if !defined(UNIX_VERSION) && !defined(WIN_VERSION) && !defined(MAC_VERSION) |
|---|
| 21 | #ifdef __MWERKS__ |
|---|
| 22 | #define MAC_VERSION |
|---|
| 23 | #elif defined __APPLE__ |
|---|
| 24 | #define MAC_VERSION |
|---|
| 25 | #else |
|---|
| 26 | #define WIN_VERSION |
|---|
| 27 | #endif |
|---|
| 28 | #endif |
|---|
| 29 | |
|---|
| 30 | /* found out that mrbayes crashes on 64 bit platform |
|---|
| 31 | especially in sumt function. If every long is substituted with |
|---|
| 32 | an int, it works. I'm going to define a SafeLong and an unsigned |
|---|
| 33 | SafeLong for 64 bit platforms... |
|---|
| 34 | Davide Cittaro - daweonline(at)gmail.com |
|---|
| 35 | */ |
|---|
| 36 | |
|---|
| 37 | //#define SSE_ENABLED |
|---|
| 38 | |
|---|
| 39 | /* This is a configuration option from the configure script. */ |
|---|
| 40 | #ifdef _64BIT |
|---|
| 41 | typedef int SafeLong; |
|---|
| 42 | #else |
|---|
| 43 | typedef long SafeLong; |
|---|
| 44 | #endif |
|---|
| 45 | |
|---|
| 46 | |
|---|
| 47 | typedef double MrBFlt; /* double used for parameter values and generally for floating point values, if set to float MPI would not work becouse of use MPI_DOUBLE*/ |
|---|
| 48 | #define MRBFLT_MAX DBL_MAX; /* maximum possible value that can be stored in MrBFlt */ |
|---|
| 49 | #define MRBFLT_MIN DBL_MIN; /* maximum possible value that can be stored in MrBFlt */ |
|---|
| 50 | #define MRBFLT_NEG_MAX (-DBL_MAX); /* maximum possible negative value that can be stored in MrBFlt */ |
|---|
| 51 | typedef float CLFlt; /* single-precision float used for cond likes (CLFlt) to increase speed and reduce memory requirement */ |
|---|
| 52 | /* set CLFlt to double if you want increased precision */ |
|---|
| 53 | /* NOTE: CLFlt = double not compatible with SSE_ENABLED */ |
|---|
| 54 | |
|---|
| 55 | |
|---|
| 56 | /* Define a compiler and vector size for the SSE code */ |
|---|
| 57 | #if defined (SSE_ENABLED) |
|---|
| 58 | #define FLOATS_PER_VEC 4 |
|---|
| 59 | #if defined (WIN_VERSION) |
|---|
| 60 | #define MS_VCPP_SSE |
|---|
| 61 | #else |
|---|
| 62 | #define GCC_SSE |
|---|
| 63 | #undef ICC_SSE |
|---|
| 64 | #endif |
|---|
| 65 | #endif |
|---|
| 66 | |
|---|
| 67 | /*#define PRINT_RATEMULTIPLIERS_CPP*/ |
|---|
| 68 | |
|---|
| 69 | |
|---|
| 70 | #if defined GCC_SSE /* gcc compiler */ |
|---|
| 71 | #define ALIGNED_MALLOC(X,Y,Z) {int tmp;tmp=posix_memalign(X,Y,Z);} /*To remove warnings*/ |
|---|
| 72 | #define ALIGNED_FREE free |
|---|
| 73 | #include <xmmintrin.h> |
|---|
| 74 | #elif defined ICC_SSE /* icc compiler */ |
|---|
| 75 | #define ALIGNED_MALLOC _mm_malloc |
|---|
| 76 | #define ALIGNED_FREE _mm_free |
|---|
| 77 | #elif defined MS_VCPP_SSE /* Visual .Net */ |
|---|
| 78 | #define ALIGNED_MALLOC _aligned_malloc |
|---|
| 79 | #define ALIGNED_FREE _aligned_free |
|---|
| 80 | #include <xmmintrin.h> |
|---|
| 81 | #else |
|---|
| 82 | #define ALIGNED_MALLOC malloc |
|---|
| 83 | #endif |
|---|
| 84 | |
|---|
| 85 | /* For comparing floating points: two values are the same if the absolute difference is less then |
|---|
| 86 | this value. |
|---|
| 87 | */ |
|---|
| 88 | #ifndef ETA |
|---|
| 89 | #define ETA (1E-30) |
|---|
| 90 | #endif |
|---|
| 91 | |
|---|
| 92 | #if defined (DEBUGOUTPUT) |
|---|
| 93 | #define DEBUG(fmt, arg) printf("%s:%d ",__FILE__,__LINE__);printf(fmt,arg); |
|---|
| 94 | #else |
|---|
| 95 | #define DEBUG(a,b) |
|---|
| 96 | #endif |
|---|
| 97 | |
|---|
| 98 | #if defined (MPI_ENABLED) |
|---|
| 99 | #include "mpi.h" |
|---|
| 100 | #endif |
|---|
| 101 | |
|---|
| 102 | #if defined (BEAGLE_ENABLED) |
|---|
| 103 | #include "libhmsbeagle/beagle.h" |
|---|
| 104 | #endif |
|---|
| 105 | |
|---|
| 106 | #define RELEASE |
|---|
| 107 | #ifdef RELEASE |
|---|
| 108 | #define VERSION_NUMBER "3.2.1" |
|---|
| 109 | #else |
|---|
| 110 | #define VERSION_NUMBER "3.2.1-svn" |
|---|
| 111 | #endif |
|---|
| 112 | |
|---|
| 113 | /* TEMPSTRSIZE determines size of temporary sprintf buffer (for SafeSprintf) */ |
|---|
| 114 | /* A patch was sent in by Allen Smith for SafeSprintf, but I could not get |
|---|
| 115 | it compiled on SGI IRIX 6.5 (too old?) with _xpg5_vsnprintf undefined. |
|---|
| 116 | The code below is a hack so SafeSprintf never has to reallocate memory. |
|---|
| 117 | \todo fix |
|---|
| 118 | */ |
|---|
| 119 | #ifdef __sgi |
|---|
| 120 | #define TEMPSTRSIZE 1000 |
|---|
| 121 | #else |
|---|
| 122 | #define TEMPSTRSIZE 200 |
|---|
| 123 | #endif |
|---|
| 124 | |
|---|
| 125 | #undef NO_ERROR |
|---|
| 126 | #undef ERROR |
|---|
| 127 | #define NO_ERROR 0 |
|---|
| 128 | #define ERROR 1 |
|---|
| 129 | #define NO_ERROR_QUIT 2 |
|---|
| 130 | #define ABORT 3 |
|---|
| 131 | #define SKIP_COMMAND 4 |
|---|
| 132 | |
|---|
| 133 | #undef FALSE |
|---|
| 134 | #undef TRUE |
|---|
| 135 | #define FALSE 0 |
|---|
| 136 | #define TRUE 1 |
|---|
| 137 | |
|---|
| 138 | #define NO 0 |
|---|
| 139 | #define YES 1 |
|---|
| 140 | |
|---|
| 141 | #define UP 0 |
|---|
| 142 | #define DOWN 1 |
|---|
| 143 | |
|---|
| 144 | #define UPPER 0 |
|---|
| 145 | #define MIDDLE 1 |
|---|
| 146 | #define LOWER 2 |
|---|
| 147 | |
|---|
| 148 | #define NONINTERACTIVE 0 |
|---|
| 149 | #define INTERACTIVE 1 |
|---|
| 150 | |
|---|
| 151 | #define STANDARD_USER 1 |
|---|
| 152 | #define DEVELOPER 3 |
|---|
| 153 | |
|---|
| 154 | #define DIFFERENT 0 |
|---|
| 155 | #define SAME 1 |
|---|
| 156 | #define CONSISTENT_WITH 2 |
|---|
| 157 | |
|---|
| 158 | #define LINETERM_UNIX 0 |
|---|
| 159 | #define LINETERM_MAC 1 |
|---|
| 160 | #define LINETERM_DOS 2 |
|---|
| 161 | |
|---|
| 162 | #define SCREENWIDTH 60 |
|---|
| 163 | #define SCREENWIDTH2 61 |
|---|
| 164 | |
|---|
| 165 | #define AVGSTDDEV 0 |
|---|
| 166 | #define MAXSTDDEV 1 |
|---|
| 167 | |
|---|
| 168 | #define NONE 0 |
|---|
| 169 | #define DNA 1 |
|---|
| 170 | #define RNA 2 |
|---|
| 171 | #define PROTEIN 3 |
|---|
| 172 | #define RESTRICTION 4 |
|---|
| 173 | #define STANDARD 5 |
|---|
| 174 | #define MIXED 6 |
|---|
| 175 | #define CONTINUOUS 7 |
|---|
| 176 | |
|---|
| 177 | #define AAMODEL_POISSON 0 |
|---|
| 178 | #define AAMODEL_JONES 1 |
|---|
| 179 | #define AAMODEL_DAY 2 |
|---|
| 180 | #define AAMODEL_MTREV 3 |
|---|
| 181 | #define AAMODEL_MTMAM 4 |
|---|
| 182 | #define AAMODEL_WAG 5 |
|---|
| 183 | #define AAMODEL_RTREV 6 |
|---|
| 184 | #define AAMODEL_CPREV 7 |
|---|
| 185 | #define AAMODEL_VT 8 |
|---|
| 186 | #define AAMODEL_BLOSUM 9 |
|---|
| 187 | #define AAMODEL_EQ 10 |
|---|
| 188 | #define AAMODEL_GTR 11 /* aa models with free parameters must be listed last */ |
|---|
| 189 | |
|---|
| 190 | #define NUCMODEL_4BY4 0 |
|---|
| 191 | #define NUCMODEL_DOUBLET 1 |
|---|
| 192 | #define NUCMODEL_CODON 2 |
|---|
| 193 | #define NUCMODEL_AA 3 |
|---|
| 194 | |
|---|
| 195 | #define NST_MIXED -1 /* anything other than 1, 2, or 6 */ |
|---|
| 196 | |
|---|
| 197 | #define MISSING 10000000 |
|---|
| 198 | #define GAP 10000001 |
|---|
| 199 | |
|---|
| 200 | #define UNORD 0 |
|---|
| 201 | #define ORD 1 |
|---|
| 202 | #define DOLLO 2 |
|---|
| 203 | #define IRREV 3 |
|---|
| 204 | |
|---|
| 205 | #define IN_CMD 0 |
|---|
| 206 | #define IN_FILE 1 |
|---|
| 207 | |
|---|
| 208 | #define NOTHING 0 |
|---|
| 209 | #define COMMAND 1 |
|---|
| 210 | #define PARAMETER 2 |
|---|
| 211 | #define EQUALSIGN 3 |
|---|
| 212 | #define COLON 4 |
|---|
| 213 | #define SEMICOLON 5 |
|---|
| 214 | #define COMMA 6 |
|---|
| 215 | #define POUNDSIGN 7 |
|---|
| 216 | #define QUESTIONMARK 8 |
|---|
| 217 | #define DASH 9 |
|---|
| 218 | #define LEFTPAR 10 |
|---|
| 219 | #define RIGHTPAR 11 |
|---|
| 220 | #define LEFTCOMMENT 12 |
|---|
| 221 | #define RIGHTCOMMENT 13 |
|---|
| 222 | #define ALPHA 14 |
|---|
| 223 | #define NUMBER 15 |
|---|
| 224 | #define RETURNSYMBOL 16 |
|---|
| 225 | #define ASTERISK 17 |
|---|
| 226 | #define BACKSLASH 18 |
|---|
| 227 | #define FORWARDSLASH 19 |
|---|
| 228 | #define EXCLAMATIONMARK 20 |
|---|
| 229 | #define PERCENT 21 |
|---|
| 230 | #define QUOTATIONMARK 22 |
|---|
| 231 | #define WEIRD 23 |
|---|
| 232 | #define UNKNOWN_TOKEN_TYPE 24 |
|---|
| 233 | #define LEFTCURL 25 |
|---|
| 234 | #define RIGHTCURL 26 |
|---|
| 235 | #define DOLLAR 27 |
|---|
| 236 | #define AMPERSAND 28 |
|---|
| 237 | |
|---|
| 238 | #define MAX_Q_RATE 100.0f |
|---|
| 239 | #define MIN_SHAPE_PARAM 0.0001f |
|---|
| 240 | #define MAX_SHAPE_PARAM 200.0f |
|---|
| 241 | #define MAX_SITE_RATE 10.0f |
|---|
| 242 | #define MAX_GAMMA_CATS 20 |
|---|
| 243 | #define MAX_GAMMA_CATS_SQUARED 400 |
|---|
| 244 | #define BRLENS_MIN 0.00000001f |
|---|
| 245 | #define BRLENS_MAX 100.0f |
|---|
| 246 | #define RELBRLENS_MIN 0.00000001f |
|---|
| 247 | #define RELBRLENS_MAX 100.0f |
|---|
| 248 | #define KAPPA_MIN 0.01f |
|---|
| 249 | #define KAPPA_MAX 1000.0f |
|---|
| 250 | #define GROWTH_MIN -1000000.0f |
|---|
| 251 | #define GROWTH_MAX 1000000.0f |
|---|
| 252 | #define RATE_MIN 0.000001f |
|---|
| 253 | #define RATE_MAX 100.0f |
|---|
| 254 | #define CPPRATEMULTIPLIER_MIN 0.001f |
|---|
| 255 | #define CPPRATEMULTIPLIER_MAX 1000.0f |
|---|
| 256 | #define SYMPI_MIN 0.000001f |
|---|
| 257 | #define SYMPI_MAX 100.0f |
|---|
| 258 | #define ALPHA_MIN 0.0001f |
|---|
| 259 | #define ALPHA_MAX 10000.0f |
|---|
| 260 | #define DIR_MIN 0.000001f |
|---|
| 261 | #define PI_MIN 0.000001f |
|---|
| 262 | #define OFFSETEXPLAMBDA_MIN 0.000001f |
|---|
| 263 | #define OFFSETEXPLAMBDA_MAX 100000.0f |
|---|
| 264 | #define TREEHEIGHT_MIN 0.00000001f |
|---|
| 265 | #define TREEHEIGHT_MAX 1000.0f |
|---|
| 266 | #define TREEAGE_MIN 0.00000001f |
|---|
| 267 | #define TREEAGE_MAX 1000000.0f |
|---|
| 268 | #define CPPLAMBDA_MIN 0.00001f |
|---|
| 269 | #define CPPLAMBDA_MAX 100.0f |
|---|
| 270 | #define TK02VAR_MIN 0.00001f |
|---|
| 271 | #define TK02VAR_MAX 100000.0f |
|---|
| 272 | #define IGRVAR_MIN 0.00001f |
|---|
| 273 | #define IGRVAR_MAX 100000.0f |
|---|
| 274 | #define OMEGA_MAX 1000000.0f |
|---|
| 275 | |
|---|
| 276 | #define POS_MIN 1E-25f |
|---|
| 277 | #define POS_MAX 1E25f |
|---|
| 278 | #define POS_INFINITY 1E25f |
|---|
| 279 | #define NEG_INFINITY -1000000.0f |
|---|
| 280 | #define POSREAL_MIN 1E-25f |
|---|
| 281 | #define POSREAL_MAX 1E25f |
|---|
| 282 | |
|---|
| 283 | #define CMD_STRING_LENGTH 100000 |
|---|
| 284 | |
|---|
| 285 | #define pos(i,j,n) ((i)*(n)+(j)) |
|---|
| 286 | |
|---|
| 287 | #define NUM_ALLOCS 92 |
|---|
| 288 | |
|---|
| 289 | #define ALLOC_MATRIX 0 |
|---|
| 290 | #define ALLOC_CHARINFO 2 |
|---|
| 291 | #define ALLOC_CHARSETS 3 |
|---|
| 292 | #define ALLOC_TAXA 4 |
|---|
| 293 | #define ALLOC_TMPSET 5 |
|---|
| 294 | #define ALLOC_PARTITIONS 6 |
|---|
| 295 | #define ALLOC_PARTITIONVARS 7 |
|---|
| 296 | #define ALLOC_TAXASETS 8 |
|---|
| 297 | #define ALLOC_CONSTRAINTS 9 |
|---|
| 298 | #define ALLOC_USERTREE 10 |
|---|
| 299 | #define ALLOC_SUMTPARAMS 11 |
|---|
| 300 | #define ALLOC_TERMSTATE 12 |
|---|
| 301 | #define ALLOC_ISPARTAMBIG 13 |
|---|
| 302 | #define ALLOC_AVAILNODES 25 |
|---|
| 303 | #define ALLOC_AVAILINDICES 26 |
|---|
| 304 | #define ALLOC_CURLNL 28 |
|---|
| 305 | #define ALLOC_CURLNPR 29 |
|---|
| 306 | #define ALLOC_CHAINID 30 |
|---|
| 307 | #define ALLOC_PARAMS 31 |
|---|
| 308 | #define ALLOC_TREE 32 |
|---|
| 309 | #define ALLOC_NODES 33 |
|---|
| 310 | #define ALLOC_LOCTAXANAMES 34 |
|---|
| 311 | #define ALLOC_COMPMATRIX 39 |
|---|
| 312 | #define ALLOC_NUMSITESOFPAT 40 |
|---|
| 313 | #define ALLOC_COMPCOLPOS 41 |
|---|
| 314 | #define ALLOC_COMPCHARPOS 42 |
|---|
| 315 | #define ALLOC_ORIGCHAR 43 |
|---|
| 316 | #define ALLOC_PARAMVALUES 46 |
|---|
| 317 | #define ALLOC_MCMCTREES 47 |
|---|
| 318 | #define ALLOC_MOVES 48 |
|---|
| 319 | #define ALLOC_PRELIKES 52 |
|---|
| 320 | #define ALLOC_SITEJUMP 54 |
|---|
| 321 | #define ALLOC_MARKOVTIS 55 |
|---|
| 322 | #define ALLOC_RATEPROBS 56 |
|---|
| 323 | #define ALLOC_STDTYPE 57 |
|---|
| 324 | #define ALLOC_PACKEDTREES 58 |
|---|
| 325 | #define ALLOC_SUMPSTRING 62 |
|---|
| 326 | #define ALLOC_SUMPINFO 63 |
|---|
| 327 | #define ALLOC_SWAPINFO 64 |
|---|
| 328 | #define ALLOC_SYMPIINDEX 65 |
|---|
| 329 | #define ALLOC_POSSELPROBS 66 |
|---|
| 330 | #define ALLOC_PBF 68 |
|---|
| 331 | #define ALLOC_LOCALTAXONCALIBRATION 69 |
|---|
| 332 | #define ALLOC_SPR_PARSSETS 72 |
|---|
| 333 | #define ALLOC_PFCOUNTERS 74 |
|---|
| 334 | #define ALLOC_FILEPOINTERS 75 |
|---|
| 335 | #define ALLOC_STATS 76 |
|---|
| 336 | #define ALLOC_DIAGNTREE 77 |
|---|
| 337 | #define ALLOC_USEDMOVES 82 |
|---|
| 338 | #define ALLOC_MODEL 83 |
|---|
| 339 | #define ALLOC_STDSTATEFREQS 84 |
|---|
| 340 | #define ALLOC_PRINTPARAM 85 |
|---|
| 341 | #define ALLOC_TREELIST 86 |
|---|
| 342 | #define ALLOC_TFILEPOS 87 |
|---|
| 343 | #define ALLOC_BEST 88 |
|---|
| 344 | #define ALLOC_SPECIESPARTITIONS 89 |
|---|
| 345 | #define ALLOC_SS 90 |
|---|
| 346 | |
|---|
| 347 | |
|---|
| 348 | #define LINKED 0 |
|---|
| 349 | #define UNLINKED 1 |
|---|
| 350 | |
|---|
| 351 | /*paramType*/ |
|---|
| 352 | #define NUM_LINKED 28 |
|---|
| 353 | #define P_TRATIO 0 |
|---|
| 354 | #define P_REVMAT 1 |
|---|
| 355 | #define P_OMEGA 2 |
|---|
| 356 | #define P_PI 3 |
|---|
| 357 | #define P_SHAPE 4 |
|---|
| 358 | #define P_PINVAR 5 |
|---|
| 359 | #define P_CORREL 6 |
|---|
| 360 | #define P_SWITCH 7 |
|---|
| 361 | #define P_RATEMULT 8 |
|---|
| 362 | #define P_TOPOLOGY 9 |
|---|
| 363 | #define P_BRLENS 10 |
|---|
| 364 | #define P_SPECRATE 11 |
|---|
| 365 | #define P_EXTRATE 12 |
|---|
| 366 | #define P_POPSIZE 13 |
|---|
| 367 | #define P_AAMODEL 14 |
|---|
| 368 | #define P_BRCORR 15 |
|---|
| 369 | #define P_BRSIGMA 16 |
|---|
| 370 | #define P_GROWTH 17 |
|---|
| 371 | #define P_CPPMULTDEV 18 |
|---|
| 372 | #define P_CPPRATE 19 |
|---|
| 373 | #define P_TK02VAR 20 |
|---|
| 374 | #define P_CPPEVENTS 21 |
|---|
| 375 | #define P_TK02BRANCHRATES 22 |
|---|
| 376 | #define P_IGRVAR 23 |
|---|
| 377 | #define P_IGRBRANCHLENS 24 |
|---|
| 378 | #define P_CLOCKRATE 25 |
|---|
| 379 | #define P_SPECIESTREE 26 |
|---|
| 380 | #define P_GENETREERATE 27 |
|---|
| 381 | |
|---|
| 382 | /* NOTE: If you add another parameter, change NUM_LINKED */ |
|---|
| 383 | |
|---|
| 384 | #define CPPm 0 /* CPP rate multipliers */ |
|---|
| 385 | #define CPPi 1 /* CPP independent rates */ |
|---|
| 386 | |
|---|
| 387 | #define MAX_NUM_USERTREES 200 /* maximum number of user trees MrBayes will read */ |
|---|
| 388 | #define MAX_CHAINS 256 /* maximum numbder of chains you can run actually only half of it becouse of m->lnLike[MAX_CHAINS] */ |
|---|
| 389 | |
|---|
| 390 | //#define PARAM_NAME_SIZE 400 |
|---|
| 391 | |
|---|
| 392 | typedef void * VoidPtr; |
|---|
| 393 | typedef int (*CmdFxn)(void); |
|---|
| 394 | typedef int (*ParmFxn)(char *, char *); |
|---|
| 395 | |
|---|
| 396 | typedef struct |
|---|
| 397 | { |
|---|
| 398 | MrBFlt sum; /* sum of standard deviations */ |
|---|
| 399 | MrBFlt max; /* maximum standard deviation */ |
|---|
| 400 | MrBFlt numPartitions; |
|---|
| 401 | MrBFlt numSamples; |
|---|
| 402 | MrBFlt avgStdDev; |
|---|
| 403 | MrBFlt **pair; |
|---|
| 404 | } STATS; |
|---|
| 405 | |
|---|
| 406 | /* enumeration for calibration prior */ |
|---|
| 407 | enum CALPRIOR |
|---|
| 408 | { |
|---|
| 409 | unconstrained, |
|---|
| 410 | fixed, |
|---|
| 411 | offsetExponential, |
|---|
| 412 | uniform |
|---|
| 413 | }; |
|---|
| 414 | |
|---|
| 415 | enum ConstraintType |
|---|
| 416 | { |
|---|
| 417 | PARTIAL, |
|---|
| 418 | NEGATIVE, |
|---|
| 419 | HARD |
|---|
| 420 | }; |
|---|
| 421 | |
|---|
| 422 | /* typedef for calibration */ |
|---|
| 423 | typedef struct calibration |
|---|
| 424 | { |
|---|
| 425 | char name[100]; |
|---|
| 426 | enum CALPRIOR prior; |
|---|
| 427 | MrBFlt max; |
|---|
| 428 | MrBFlt min; |
|---|
| 429 | MrBFlt offset; |
|---|
| 430 | MrBFlt lambda; |
|---|
| 431 | MrBFlt age; |
|---|
| 432 | } |
|---|
| 433 | Calibration; |
|---|
| 434 | |
|---|
| 435 | /* typedef for tree (topology) list element */ |
|---|
| 436 | typedef struct element { |
|---|
| 437 | struct element *next; |
|---|
| 438 | int *order; |
|---|
| 439 | } TreeListElement; |
|---|
| 440 | |
|---|
| 441 | /* typedef for list of trees (topologies) */ |
|---|
| 442 | typedef struct { |
|---|
| 443 | TreeListElement *first; |
|---|
| 444 | TreeListElement *last; |
|---|
| 445 | } TreeList; |
|---|
| 446 | |
|---|
| 447 | /* typedef for packed tree */ |
|---|
| 448 | typedef struct { |
|---|
| 449 | int *order; |
|---|
| 450 | MrBFlt *brlens; |
|---|
| 451 | } PackedTree; |
|---|
| 452 | |
|---|
| 453 | /* typedef for binary tree node */ |
|---|
| 454 | /* NOTE: Any variable added here must also be copied in CopyTrees */ |
|---|
| 455 | typedef struct node |
|---|
| 456 | { |
|---|
| 457 | char *label; /*!< name of node if tip */ |
|---|
| 458 | struct node *left, *right, *anc; /*!< pointers to adjacent nodes */ |
|---|
| 459 | int memoryIndex; /*!< memory index (do not change) */ |
|---|
| 460 | int index; /*!< index to node (0 to numLocalTaxa for tips) */ |
|---|
| 461 | int upDateCl; /*!< cond likes need update? */ |
|---|
| 462 | int upDateTi; /*!< transition probs need update? */ |
|---|
| 463 | int marked, x, y; /*!< scratch variables */ |
|---|
| 464 | int scalerNode; /*!< is node scaling cond likes? */ |
|---|
| 465 | int isLocked; /*!< is node locked? */ |
|---|
| 466 | int lockID; /*!< id of lock */ |
|---|
| 467 | int isDated; /*!< is node dated (calibrated)? */ |
|---|
| 468 | SafeLong *partition; /*!< pointer to bitfield describing splits */ |
|---|
| 469 | MrBFlt length; /*!< length of pending branch */ |
|---|
| 470 | MrBFlt nodeDepth; /*!< node depth (height) */ |
|---|
| 471 | MrBFlt age; /*!< age of node */ |
|---|
| 472 | MrBFlt d; /*!< scratch variable */ |
|---|
| 473 | Calibration *calibration; /*!< pointer to calibration data */ |
|---|
| 474 | } |
|---|
| 475 | TreeNode; |
|---|
| 476 | |
|---|
| 477 | |
|---|
| 478 | /* typedef for binary tree */ |
|---|
| 479 | typedef struct |
|---|
| 480 | { |
|---|
| 481 | char name[100]; /*!< name of tree */ |
|---|
| 482 | int memNodes; /*!< number of allocated nodes (do not exceed!) */ |
|---|
| 483 | int nNodes; /*!< number of nodes in tree (including lower root in rooted trees) */ |
|---|
| 484 | int nIntNodes; /*!< number of interior nodes in tree (excluding lower root in rooted trees) */ |
|---|
| 485 | int isRooted; /*!< is tree rooted? */ |
|---|
| 486 | int isClock; /*!< is tree clock? */ |
|---|
| 487 | int isCalibrated; /*!< is tree calibrated? */ |
|---|
| 488 | int nRelParts; /*!< number of relevant partitions */ |
|---|
| 489 | int *relParts; /*!< pointer to relevant partitions */ |
|---|
| 490 | int checkConstraints; /*!< does tree have constraints? */ |
|---|
| 491 | int nConstraints; /*!< number of constraints */ |
|---|
| 492 | int *constraints; /*!< pointer to constraints */ |
|---|
| 493 | int nLocks; /*!< number of constrained (locked) nodes */ |
|---|
| 494 | TreeNode **allDownPass; /*!< downpass array of all nodes */ |
|---|
| 495 | TreeNode **intDownPass; /*!< downpass array of interior nodes (including upper but excluding lower root in rooted trees) */ |
|---|
| 496 | TreeNode *root; /*!< pointer to root (lower root in rooted trees) */ |
|---|
| 497 | TreeNode *nodes; /*!< array containing the nodes */ |
|---|
| 498 | SafeLong *bitsets; /*!< pointer to bitsets describing splits */ |
|---|
| 499 | SafeLong *flags; /*!< pointer to cond like flags */ |
|---|
| 500 | 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 */ |
|---|
| 501 | } |
|---|
| 502 | Tree; |
|---|
| 503 | |
|---|
| 504 | |
|---|
| 505 | /* typedef for node in polytomous tree */ |
|---|
| 506 | typedef struct pNode |
|---|
| 507 | { |
|---|
| 508 | char label[100]; /*!< name of node if terminal */ |
|---|
| 509 | struct pNode *left, *sib, *anc; /*!< pointers to adjacent nodes */ |
|---|
| 510 | int x, y, mark; /*!< scratch variables */ |
|---|
| 511 | int partitionIndex; /*!< partition index in sumt (scratch) */ |
|---|
| 512 | int index; /*!< index of node (if < numLocalTaxa = local taxon index) */ |
|---|
| 513 | int memoryIndex; /*!< immutable index of memory position */ |
|---|
| 514 | int isLocked; /*!< is the node locked? */ |
|---|
| 515 | int lockID; /*!< id of lock */ |
|---|
| 516 | int isDated; /*!< is node dated? */ |
|---|
| 517 | MrBFlt length; /*!< age of node */ |
|---|
| 518 | MrBFlt depth; /*!< depth (height) of node */ |
|---|
| 519 | MrBFlt age; /*!< age of node */ |
|---|
| 520 | MrBFlt support, f; /*!< scratch variables */ |
|---|
| 521 | SafeLong *partition; /*!< pointer to partition (split) bitset */ |
|---|
| 522 | Calibration *calibration; /*!< pointer to dating of node */ |
|---|
| 523 | } |
|---|
| 524 | PolyNode; |
|---|
| 525 | |
|---|
| 526 | |
|---|
| 527 | /* typedef for polytomous tree */ |
|---|
| 528 | typedef struct |
|---|
| 529 | { |
|---|
| 530 | char name[100]; /*!< name of tree */ |
|---|
| 531 | int memNodes; /*!< number of allocated nodes; do not exceed! */ |
|---|
| 532 | int nNodes; /*!< number of nodes in tree */ |
|---|
| 533 | int nIntNodes; /*!< number of interior nodes in tree */ |
|---|
| 534 | PolyNode **allDownPass; /*!< downpass array over all nodes */ |
|---|
| 535 | PolyNode **intDownPass; /*!< downpass array over interior nodes */ |
|---|
| 536 | PolyNode *root; /*!< pointer to root (lower for rooted trees */ |
|---|
| 537 | PolyNode *nodes; /*!< array holding the tree nodes */ |
|---|
| 538 | SafeLong *bitsets; /*!< bits describing partitions (splits) */ |
|---|
| 539 | int nBSets; /*!< number of effective branch length sets */ |
|---|
| 540 | int nESets; /*!< number of breakpoint rate sets */ |
|---|
| 541 | char **bSetName; /*!< names of effective branch length sets */ |
|---|
| 542 | char **eSetName; /*!< names of breakpoint rate sets */ |
|---|
| 543 | int **nEvents; /*!< number of branch events of bp rate set */ |
|---|
| 544 | MrBFlt ***position; /*!< position of branch events */ |
|---|
| 545 | MrBFlt ***rateMult; /*!< parameter of branch events */ |
|---|
| 546 | MrBFlt **effectiveBrLen; /*!< effective branch lengths of ebl set */ |
|---|
| 547 | int brlensDef; /*!< are brlens defined ? */ |
|---|
| 548 | int isRooted; /*!< is tree rooted? */ |
|---|
| 549 | int isClock; /*!< is tree clock? */ |
|---|
| 550 | int isCalibrated; /*!< is tree calibrated? */ |
|---|
| 551 | int isRelaxed; /*!< is tree relaxed? */ |
|---|
| 552 | MrBFlt clockRate; /*!< clock rate */ |
|---|
| 553 | int popSizeSet; /*!< does tree have a population size set? */ |
|---|
| 554 | MrBFlt *popSize; /*!< the population size */ |
|---|
| 555 | char *popSizeSetName; /*!< name of the population size set */ |
|---|
| 556 | } |
|---|
| 557 | PolyTree; |
|---|
| 558 | |
|---|
| 559 | /* typedef for a ln prior prob fxn */ |
|---|
| 560 | typedef MrBFlt (*LnPriorProbFxn)(MrBFlt val, MrBFlt *priorParams); |
|---|
| 561 | |
|---|
| 562 | /* typedef for a ln prior prob ratio fxn */ |
|---|
| 563 | typedef MrBFlt (*LnPriorRatioFxn)(MrBFlt newVal, MrBFlt oldVal, MrBFlt *priorParams); |
|---|
| 564 | |
|---|
| 565 | /* struct for holding model parameter info for the mcmc run */ |
|---|
| 566 | typedef struct param |
|---|
| 567 | { |
|---|
| 568 | int index; /* index to the parameter (0, 1, 2, ...) */ |
|---|
| 569 | int paramType; /* the type of the parameter */ |
|---|
| 570 | int paramId; /* unique ID for parameter x prior combination */ |
|---|
| 571 | MrBFlt *values; /* main values of parameter */ |
|---|
| 572 | MrBFlt *subValues; /* subvalues of parameter */ |
|---|
| 573 | int *intValues; /* integer values (model index/growth fxn) */ |
|---|
| 574 | int nValues; /* number of values */ |
|---|
| 575 | int nSubValues; /* number of subvalues */ |
|---|
| 576 | int nIntValues; /* number of intvalues */ |
|---|
| 577 | MrBFlt min; /* minimum value of parameter */ |
|---|
| 578 | MrBFlt max; /* maximum value of parameter */ |
|---|
| 579 | int *relParts; /* pointer to relevant divisions */ |
|---|
| 580 | int nRelParts; /* number of relevant divisions */ |
|---|
| 581 | int upDate; /* update flag (for copying) */ |
|---|
| 582 | struct param **subParams; /* pointers to subparams (for topology) */ |
|---|
| 583 | int nSubParams; /* number of subparams */ |
|---|
| 584 | Tree **tree; /* pointer to tree ptrs (for brlens & topology) */ |
|---|
| 585 | int treeIndex; /* index to first tree in mcmcTree */ |
|---|
| 586 | int hasBinaryStd; /* has binary standard chars */ |
|---|
| 587 | int *sympiBsIndex; /* pointer to sympi bsIndex (std chars) */ |
|---|
| 588 | int *sympinStates; /* pointer to sympi nStates (std chars) */ |
|---|
| 589 | int *sympiCType; /* pointer to sympi cType (std chars) */ |
|---|
| 590 | int nSympi; /* number of sympis */ |
|---|
| 591 | int printParam; /* whether parameter should be printed */ |
|---|
| 592 | int nPrintSubParams; /* number of subparams that should be printed */ |
|---|
| 593 | char *paramHeader; /* a string holding header for param values */ |
|---|
| 594 | char *name; /* string holding name of parameter */ |
|---|
| 595 | char *paramTypeName; /* pointer to description of parameter type */ |
|---|
| 596 | int checkConstraints; /* is tree parameter constrained? */ |
|---|
| 597 | int fill; /* flags whether the parameter should be filled */ |
|---|
| 598 | int nStdStateFreqs; /* number of std state frequencies */ |
|---|
| 599 | MrBFlt *stdStateFreqs; /* pointer to std state frequencies */ |
|---|
| 600 | int **nEvents; /* number of branch events for Cpp model */ |
|---|
| 601 | /* nEvents[0..2*numCains][0..numNodes=2*numTaxa]*/ |
|---|
| 602 | MrBFlt ***position; /* event positions for Cpp relaxed clock model */ |
|---|
| 603 | MrBFlt ***rateMult; /* rate multipliers for Cpp relaxed clock model */ |
|---|
| 604 | int affectsLikelihood; /* does parameter directly influence likelihood? */ |
|---|
| 605 | MrBFlt* priorParams; /* pointer to the prior parameters */ |
|---|
| 606 | LnPriorProbFxn LnPriorProb; /* ln prior prob function */ |
|---|
| 607 | LnPriorRatioFxn LnPriorRatio; /* ln prior prob ratio function */ |
|---|
| 608 | } Param; |
|---|
| 609 | |
|---|
| 610 | |
|---|
| 611 | #if defined(THREADS_ENABLED) |
|---|
| 612 | #include <pthread.h> |
|---|
| 613 | |
|---|
| 614 | typedef struct s_launch_struct |
|---|
| 615 | { |
|---|
| 616 | int chain; |
|---|
| 617 | int division; |
|---|
| 618 | MrBFlt* lnL; |
|---|
| 619 | } LaunchStruct; |
|---|
| 620 | #endif |
|---|
| 621 | |
|---|
| 622 | /* parameter ID values */ |
|---|
| 623 | /* identifies unique model parameter x prior combinations */ |
|---|
| 624 | #define TRATIO_DIR 1 |
|---|
| 625 | #define TRATIO_FIX 2 |
|---|
| 626 | #define REVMAT_DIR 3 |
|---|
| 627 | #define REVMAT_FIX 4 |
|---|
| 628 | #define OMEGA_DIR 5 |
|---|
| 629 | #define OMEGA_FIX 6 |
|---|
| 630 | #define SYMPI_UNI 7 |
|---|
| 631 | #define SYMPI_UNI_MS 8 |
|---|
| 632 | #define SYMPI_EXP 9 |
|---|
| 633 | #define SYMPI_EXP_MS 10 |
|---|
| 634 | #define SYMPI_FIX 11 |
|---|
| 635 | #define SYMPI_FIX_MS 12 |
|---|
| 636 | #define SYMPI_EQUAL 13 |
|---|
| 637 | #define PI_DIR 14 |
|---|
| 638 | #define PI_USER 15 |
|---|
| 639 | #define PI_EMPIRICAL 16 |
|---|
| 640 | #define PI_EQUAL 17 |
|---|
| 641 | #define PI_FIXED 18 |
|---|
| 642 | #define SHAPE_UNI 19 |
|---|
| 643 | #define SHAPE_EXP 20 |
|---|
| 644 | #define SHAPE_FIX 21 |
|---|
| 645 | #define PINVAR_UNI 22 |
|---|
| 646 | #define PINVAR_FIX 23 |
|---|
| 647 | #define CORREL_UNI 24 |
|---|
| 648 | #define CORREL_FIX 25 |
|---|
| 649 | #define SWITCH_UNI 26 |
|---|
| 650 | #define SWITCH_EXP 27 |
|---|
| 651 | #define SWITCH_FIX 28 |
|---|
| 652 | #define RATEMULT_DIR 29 |
|---|
| 653 | #define RATEMULT_FIX 30 |
|---|
| 654 | #define TOPOLOGY_NCL_UNIFORM 31 |
|---|
| 655 | #define TOPOLOGY_NCL_CONSTRAINED 32 |
|---|
| 656 | #define TOPOLOGY_NCL_FIXED 33 |
|---|
| 657 | #define TOPOLOGY_NCL_UNIFORM_HOMO 34 |
|---|
| 658 | #define TOPOLOGY_NCL_UNIFORM_HETERO 35 |
|---|
| 659 | #define TOPOLOGY_NCL_CONSTRAINED_HOMO 36 |
|---|
| 660 | #define TOPOLOGY_NCL_CONSTRAINED_HETERO 37 |
|---|
| 661 | #define TOPOLOGY_NCL_FIXED_HOMO 38 |
|---|
| 662 | #define TOPOLOGY_NCL_FIXED_HETERO 39 |
|---|
| 663 | #define TOPOLOGY_CL_UNIFORM 40 |
|---|
| 664 | #define TOPOLOGY_CL_CONSTRAINED 41 |
|---|
| 665 | #define TOPOLOGY_CL_FIXED 42 |
|---|
| 666 | #define TOPOLOGY_CCL_UNIFORM 43 |
|---|
| 667 | #define TOPOLOGY_CCL_CONSTRAINED 44 |
|---|
| 668 | #define TOPOLOGY_CCL_FIXED 45 |
|---|
| 669 | #define TOPOLOGY_PARSIMONY_UNIFORM 46 |
|---|
| 670 | #define TOPOLOGY_PARSIMONY_CONSTRAINED 47 |
|---|
| 671 | #define TOPOLOGY_PARSIMONY_FIXED 48 |
|---|
| 672 | #define BRLENS_UNI 49 |
|---|
| 673 | #define BRLENS_EXP 50 |
|---|
| 674 | #define BRLENS_FIXED 51 |
|---|
| 675 | #define BRLENS_CLOCK_UNI 52 |
|---|
| 676 | #define BRLENS_CLOCK_COAL 53 |
|---|
| 677 | #define BRLENS_CLOCK_BD 54 |
|---|
| 678 | #define BRLENS_CLOCK_FIXED 55 |
|---|
| 679 | #define BRLENS_CLOCK_SPCOAL 56 |
|---|
| 680 | #define BRLENS_PARSIMONY 57 |
|---|
| 681 | #define SPECRATE_UNI 60 |
|---|
| 682 | #define SPECRATE_EXP 61 |
|---|
| 683 | #define SPECRATE_FIX 62 |
|---|
| 684 | #define EXTRATE_BETA 63 |
|---|
| 685 | #define EXTRATE_FIX 65 |
|---|
| 686 | #define POPSIZE_UNI 66 |
|---|
| 687 | #define POPSIZE_GAMMA 67 |
|---|
| 688 | #define POPSIZE_FIX 68 |
|---|
| 689 | #define POPSIZE_NORMAL 69 |
|---|
| 690 | #define POPSIZE_LOGNORMAL 70 |
|---|
| 691 | #define AAMODEL_FIX 71 |
|---|
| 692 | #define AAMODEL_MIX 72 |
|---|
| 693 | #define GROWTH_UNI 73 |
|---|
| 694 | #define GROWTH_EXP 74 |
|---|
| 695 | #define GROWTH_FIX 75 |
|---|
| 696 | #define GROWTH_NORMAL 76 |
|---|
| 697 | #define OMEGA_BUD 77 |
|---|
| 698 | #define OMEGA_BUF 78 |
|---|
| 699 | #define OMEGA_BED 79 |
|---|
| 700 | #define OMEGA_BEF 80 |
|---|
| 701 | #define OMEGA_BFD 81 |
|---|
| 702 | #define OMEGA_BFF 82 |
|---|
| 703 | #define OMEGA_FUD 83 |
|---|
| 704 | #define OMEGA_FUF 84 |
|---|
| 705 | #define OMEGA_FED 85 |
|---|
| 706 | #define OMEGA_FEF 86 |
|---|
| 707 | #define OMEGA_FFD 87 |
|---|
| 708 | #define OMEGA_FFF 88 |
|---|
| 709 | #define OMEGA_ED 89 |
|---|
| 710 | #define OMEGA_EF 90 |
|---|
| 711 | #define OMEGA_FD 91 |
|---|
| 712 | #define OMEGA_FF 92 |
|---|
| 713 | #define OMEGA_10UUB 93 |
|---|
| 714 | #define OMEGA_10UUF 94 |
|---|
| 715 | #define OMEGA_10UEB 95 |
|---|
| 716 | #define OMEGA_10UEF 96 |
|---|
| 717 | #define OMEGA_10UFB 97 |
|---|
| 718 | #define OMEGA_10UFF 98 |
|---|
| 719 | #define OMEGA_10EUB 99 |
|---|
| 720 | #define OMEGA_10EUF 100 |
|---|
| 721 | #define OMEGA_10EEB 101 |
|---|
| 722 | #define OMEGA_10EEF 102 |
|---|
| 723 | #define OMEGA_10EFB 103 |
|---|
| 724 | #define OMEGA_10EFF 104 |
|---|
| 725 | #define OMEGA_10FUB 105 |
|---|
| 726 | #define OMEGA_10FUF 106 |
|---|
| 727 | #define OMEGA_10FEB 107 |
|---|
| 728 | #define OMEGA_10FEF 108 |
|---|
| 729 | #define OMEGA_10FFB 109 |
|---|
| 730 | #define OMEGA_10FFF 110 |
|---|
| 731 | #define CPPRATE_FIX 111 |
|---|
| 732 | #define CPPRATE_EXP 112 |
|---|
| 733 | #define CPPMULTDEV_FIX 113 |
|---|
| 734 | #define TK02VAR_FIX 114 |
|---|
| 735 | #define TK02VAR_EXP 115 |
|---|
| 736 | #define TK02VAR_UNI 116 |
|---|
| 737 | #define TOPOLOGY_RCL_UNIFORM 117 |
|---|
| 738 | #define TOPOLOGY_RCL_CONSTRAINED 118 |
|---|
| 739 | #define TOPOLOGY_RCL_FIXED 119 |
|---|
| 740 | #define TOPOLOGY_RCCL_UNIFORM 121 |
|---|
| 741 | #define TOPOLOGY_RCCL_CONSTRAINED 122 |
|---|
| 742 | #define TOPOLOGY_RCCL_FIXED 123 |
|---|
| 743 | #define TOPOLOGY_SPECIESTREE 124 |
|---|
| 744 | #define CPPEVENTS 125 |
|---|
| 745 | #define TK02BRANCHRATES 126 |
|---|
| 746 | #define TOPOLOGY_FIXED 127 |
|---|
| 747 | #define IGRVAR_FIX 129 |
|---|
| 748 | #define IGRVAR_EXP 130 |
|---|
| 749 | #define IGRVAR_UNI 131 |
|---|
| 750 | #define IGRBRANCHLENS 132 |
|---|
| 751 | #define CLOCKRATE_FIX 133 |
|---|
| 752 | #define CLOCKRATE_NORMAL 134 |
|---|
| 753 | #define CLOCKRATE_LOGNORMAL 135 |
|---|
| 754 | #define CLOCKRATE_GAMMA 136 |
|---|
| 755 | #define CLOCKRATE_EXP 137 |
|---|
| 756 | #define SPECIESTREE_UNIFORM 138 |
|---|
| 757 | #define GENETREERATEMULT_DIR 139 |
|---|
| 758 | #define REVMAT_MIX 140 |
|---|
| 759 | |
|---|
| 760 | |
|---|
| 761 | |
|---|
| 762 | #if defined (BEAGLE_ENABLED) |
|---|
| 763 | #define MB_BEAGLE_SCALE_ALWAYS 0 |
|---|
| 764 | #define MB_BEAGLE_SCALE_DYNAMIC 1 |
|---|
| 765 | #if defined (_DEBUG) |
|---|
| 766 | #define MB_PRINT_DYNAMIC_RESCALE_FAIL_STAT |
|---|
| 767 | #endif |
|---|
| 768 | #endif |
|---|
| 769 | |
|---|
| 770 | |
|---|
| 771 | /* typedef for a MoveFxn */ |
|---|
| 772 | typedef int (MoveFxn)(Param *param, int chain, SafeLong *seed, MrBFlt *lnPriorRatio, MrBFlt *lnProposalRatio, MrBFlt *mvp); |
|---|
| 773 | |
|---|
| 774 | /* typedef for an ApplicFxn */ |
|---|
| 775 | typedef int (ApplicFxn)(Param *param); |
|---|
| 776 | |
|---|
| 777 | /* typedef for an AutotuneFxn */ |
|---|
| 778 | typedef void (AutotuneFxn)(MrBFlt acceptanceRate, MrBFlt targetRate, int batch, MrBFlt *tuningParameter, MrBFlt minTuning, MrBFlt maxTuning); |
|---|
| 779 | |
|---|
| 780 | /* struct holding info on each move type that the program handles */ |
|---|
| 781 | typedef struct |
|---|
| 782 | { |
|---|
| 783 | MoveFxn *moveFxn; /* pointer to the move function */ |
|---|
| 784 | ApplicFxn *isApplicable; /* pointer to function determining whether move is applicable to a parameter */ |
|---|
| 785 | int nApplicable; /* number of relevant params */ |
|---|
| 786 | int applicableTo[40]; /* pointer to ID of relevant params */ |
|---|
| 787 | char *name; /* name of the move type */ |
|---|
| 788 | char *shortName; /* abbreviated name of the move type */ |
|---|
| 789 | char *paramName; /* name of subparameter if complex parameter */ |
|---|
| 790 | int subParams; /* are we changing subparams (brlens of topol.) */ |
|---|
| 791 | char *tuningName[3]; /* name of tuning params */ |
|---|
| 792 | char *shortTuningName[3];/* short name of tuning params */ |
|---|
| 793 | MrBFlt relProposalProb; /* default relative proposal probability */ |
|---|
| 794 | int numTuningParams; /* number of tuning parameters */ |
|---|
| 795 | MrBFlt tuningParam[3]; /* default tuning parameters for the proposal */ |
|---|
| 796 | MrBFlt minimum[3]; /* minimum values for tuning params */ |
|---|
| 797 | MrBFlt maximum[3]; /* maximum values for tuning params */ |
|---|
| 798 | int parsimonyBased; /* this move is based on parsimony (YES/NO) */ |
|---|
| 799 | int level; /* user level of this move */ |
|---|
| 800 | AutotuneFxn *Autotune; /* pointer to the autotuning function */ |
|---|
| 801 | MrBFlt targetRate; /* default target acceptance rate for autotuning */ |
|---|
| 802 | } MoveType; |
|---|
| 803 | |
|---|
| 804 | |
|---|
| 805 | /* max number of move types */ |
|---|
| 806 | #define NUM_MOVE_TYPES 100 |
|---|
| 807 | |
|---|
| 808 | |
|---|
| 809 | /* struct holding info on each move */ |
|---|
| 810 | /* Note: This allows several different moves to affect the same parameter */ |
|---|
| 811 | /* It also allows the same move to affect different parameters as before */ |
|---|
| 812 | /* This is also a good place to keep the proposal probs */ |
|---|
| 813 | typedef struct |
|---|
| 814 | { |
|---|
| 815 | char *name; /* name of the move */ |
|---|
| 816 | MoveType *moveType; /* pointer to the move type */ |
|---|
| 817 | MoveFxn *moveFxn; /* pointer to the move function */ |
|---|
| 818 | Param *parm; /* ptr to parameter the move applies to */ |
|---|
| 819 | MrBFlt *relProposalProb; /* the relative proposal probability */ |
|---|
| 820 | MrBFlt *cumProposalProb; /* the cumulative proposal probability */ |
|---|
| 821 | int *nAccepted; /* number of accepted moves */ |
|---|
| 822 | int *nTried; /* number of tried moves */ |
|---|
| 823 | int *nBatches; /* counter for autotuning rounds */ |
|---|
| 824 | int *nTotAccepted; /* total number of accepted moves */ |
|---|
| 825 | int *nTotTried; /* total number of tried moves */ |
|---|
| 826 | MrBFlt *targetRate; /* target acceptance rate for autotuning */ |
|---|
| 827 | MrBFlt *lastAcceptanceRate;/* acceptance rate in last complete batch */ |
|---|
| 828 | MrBFlt **tuningParam; /* tuning parameters for the move */ |
|---|
| 829 | } MCMCMove; |
|---|
| 830 | |
|---|
| 831 | typedef int (*LikeDownFxn)(TreeNode *, int, int); |
|---|
| 832 | typedef int (*LikeRootFxn)(TreeNode *, int, int); |
|---|
| 833 | typedef int (*LikeScalerFxn)(TreeNode *, int, int); |
|---|
| 834 | typedef int (*LikeFxn)(TreeNode *, int, int, MrBFlt *, int); |
|---|
| 835 | typedef int (*TiProbFxn)(TreeNode *, int, int); |
|---|
| 836 | typedef int (*LikeUpFxn)(TreeNode *, int, int); |
|---|
| 837 | typedef int (*PrintAncStFxn)(TreeNode *, int, int); |
|---|
| 838 | typedef int (*StateCodeFxn) (int); |
|---|
| 839 | typedef int (*PrintSiteRateFxn) (TreeNode *, int, int); |
|---|
| 840 | |
|---|
| 841 | |
|---|
| 842 | typedef struct cmdtyp |
|---|
| 843 | { |
|---|
| 844 | int cmdNumber; |
|---|
| 845 | char *string; |
|---|
| 846 | int specialCmd; |
|---|
| 847 | CmdFxn cmdFxnPtr; |
|---|
| 848 | short numParms; |
|---|
| 849 | short parmList[50]; |
|---|
| 850 | int expect; |
|---|
| 851 | char *cmdDescription; |
|---|
| 852 | int cmdUse; |
|---|
| 853 | int hiding; |
|---|
| 854 | } |
|---|
| 855 | CmdType; |
|---|
| 856 | |
|---|
| 857 | typedef struct parm |
|---|
| 858 | { |
|---|
| 859 | char *string; /* parameter name */ |
|---|
| 860 | char *valueList; /* list of values that could be input */ |
|---|
| 861 | ParmFxn fp; /* function pointer */ |
|---|
| 862 | } |
|---|
| 863 | ParmInfo, *ParmInfoPtr; |
|---|
| 864 | |
|---|
| 865 | typedef struct model |
|---|
| 866 | { |
|---|
| 867 | int dataType; /* data type for partition */ |
|---|
| 868 | int nStates; /* number of states for this type of data */ |
|---|
| 869 | int codon[64]; /* gives protein ID for each codon */ |
|---|
| 870 | int codonNucs[64][3]; /* gives triplet for each codon */ |
|---|
| 871 | int codonAAs[64]; /* gives protein ID for implemented code */ |
|---|
| 872 | |
|---|
| 873 | char nucModel[100]; /* nucleotide model used */ |
|---|
| 874 | char nst[100]; /* number of substitution types */ |
|---|
| 875 | char parsModel[100]; /* use the (so-called) parsimony model */ |
|---|
| 876 | char geneticCode[100]; /* genetic code used */ |
|---|
| 877 | char coding[100]; /* type of patterns encoded */ |
|---|
| 878 | char ploidy[100]; /* ploidy level */ |
|---|
| 879 | char omegaVar[100]; /* type of omega variation model */ |
|---|
| 880 | char ratesModel[100]; /* rates across sites model */ |
|---|
| 881 | int numGammaCats; /* number of categories for gamma approximation */ |
|---|
| 882 | char useGibbs[100]; /* flags whether Gibbs sampling of discrete gamma is used */ |
|---|
| 883 | int gibbsFreq; /* frequency of Gibbs resampling of discrete gamma */ |
|---|
| 884 | |
|---|
| 885 | int numBetaCats; /* number of categories for beta approximation */ |
|---|
| 886 | int numM10GammaCats; /* number of cats for gamma approx (M10 model) */ |
|---|
| 887 | int numM10BetaCats; /* number of cats for beta approx (M10 model) */ |
|---|
| 888 | char covarionModel[100];/* use covarion model? (yes/no) */ |
|---|
| 889 | char augmentData[100]; /* should data be augmented */ |
|---|
| 890 | |
|---|
| 891 | char tRatioPr[100]; /* prior for ti/tv rate ratio */ |
|---|
| 892 | MrBFlt tRatioFix; |
|---|
| 893 | MrBFlt tRatioDir[2]; |
|---|
| 894 | char revMatPr[100]; /* prior for GTR model */ |
|---|
| 895 | MrBFlt revMatFix[6]; |
|---|
| 896 | MrBFlt revMatDir[6]; |
|---|
| 897 | MrBFlt revMatSymDir; /* prior for mixed GTR subspace model */ |
|---|
| 898 | char aaModelPr[100]; /* prior for amino acid model */ |
|---|
| 899 | char aaModel[100]; |
|---|
| 900 | MrBFlt aaModelPrProbs[10]; |
|---|
| 901 | char aaRevMatPr[100]; /* prior for aa GTR model */ |
|---|
| 902 | MrBFlt aaRevMatFix[190]; |
|---|
| 903 | MrBFlt aaRevMatDir[190]; |
|---|
| 904 | char omegaPr[100]; /* prior for omega */ |
|---|
| 905 | MrBFlt omegaFix; |
|---|
| 906 | MrBFlt omegaDir[2]; |
|---|
| 907 | char ny98omega1pr[100]; /* prior for class 1 omega (Ny98 model) */ |
|---|
| 908 | MrBFlt ny98omega1Fixed; |
|---|
| 909 | MrBFlt ny98omega1Beta[2]; |
|---|
| 910 | char ny98omega3pr[100]; /* prior for class 3 omega (Ny98 model) */ |
|---|
| 911 | MrBFlt ny98omega3Fixed; |
|---|
| 912 | MrBFlt ny98omega3Uni[2]; |
|---|
| 913 | MrBFlt ny98omega3Exp; |
|---|
| 914 | char m3omegapr[100]; /* prior for all three omegas (M3 model) */ |
|---|
| 915 | MrBFlt m3omegaFixed[3]; |
|---|
| 916 | char m10betapr[100]; /* prior for omega variation (M10 model) */ |
|---|
| 917 | char m10gammapr[100]; |
|---|
| 918 | MrBFlt m10betaExp; |
|---|
| 919 | MrBFlt m10betaUni[2]; |
|---|
| 920 | MrBFlt m10betaFix[2]; |
|---|
| 921 | MrBFlt m10gammaExp; |
|---|
| 922 | MrBFlt m10gammaUni[2]; |
|---|
| 923 | MrBFlt m10gammaFix[2]; |
|---|
| 924 | char codonCatFreqPr[100];/* prior for selection cat frequencies */ |
|---|
| 925 | MrBFlt codonCatFreqFix[3]; |
|---|
| 926 | MrBFlt codonCatDir[3]; |
|---|
| 927 | char stateFreqPr[100]; /* prior for character state frequencies */ |
|---|
| 928 | MrBFlt stateFreqsFix[200]; |
|---|
| 929 | MrBFlt stateFreqsDir[200]; |
|---|
| 930 | char stateFreqsFixType[100]; |
|---|
| 931 | int numDirParams; |
|---|
| 932 | char shapePr[100]; /* prior for gamma shape parameter */ |
|---|
| 933 | MrBFlt shapeFix; |
|---|
| 934 | MrBFlt shapeUni[2]; |
|---|
| 935 | MrBFlt shapeExp; |
|---|
| 936 | char pInvarPr[100]; /* prior for proportion of invariable sites */ |
|---|
| 937 | MrBFlt pInvarFix; |
|---|
| 938 | MrBFlt pInvarUni[2]; |
|---|
| 939 | char adGammaCorPr[100]; /* prior for correlation param of adGamma model */ |
|---|
| 940 | MrBFlt corrFix; |
|---|
| 941 | MrBFlt corrUni[2]; |
|---|
| 942 | char covSwitchPr[100]; /* prior for switching rates of covarion model */ |
|---|
| 943 | MrBFlt covswitchFix[2]; |
|---|
| 944 | MrBFlt covswitchUni[2]; |
|---|
| 945 | MrBFlt covswitchExp; |
|---|
| 946 | char symPiPr[100]; /* prior for pi when unidentifiable states used */ |
|---|
| 947 | MrBFlt symBetaFix; |
|---|
| 948 | MrBFlt symBetaUni[2]; |
|---|
| 949 | MrBFlt symBetaExp; |
|---|
| 950 | char ratePr[100]; /* prior on rate for a partition */ |
|---|
| 951 | MrBFlt ratePrDir; |
|---|
| 952 | char brownCorPr[100]; /* prior for correlation of Brownian model */ |
|---|
| 953 | MrBFlt brownCorrFix; |
|---|
| 954 | MrBFlt brownCorrUni[2]; |
|---|
| 955 | char brownScalesPr[100];/* prior for scales of Brownian model */ |
|---|
| 956 | MrBFlt brownScalesFix; |
|---|
| 957 | MrBFlt brownScalesUni[2]; |
|---|
| 958 | MrBFlt brownScalesGamma[2]; |
|---|
| 959 | MrBFlt brownScalesGammaMean; |
|---|
| 960 | |
|---|
| 961 | char topologyPr[100]; /* prior for tree topology */ |
|---|
| 962 | int topologyFix; /* user tree index for fixed topology */ |
|---|
| 963 | int *activeConstraints;/* which constraints are active? */ |
|---|
| 964 | int numActiveConstraints; |
|---|
| 965 | int numActiveLocks; |
|---|
| 966 | char brlensPr[100]; /* prior on branch lengths */ |
|---|
| 967 | int brlensFix; /* user tree index for fixed brlens */ |
|---|
| 968 | MrBFlt brlensUni[2]; |
|---|
| 969 | MrBFlt brlensExp; |
|---|
| 970 | char speciesTreeBrlensPr[100]; /* prior on branch lengths of species tree */ |
|---|
| 971 | char unconstrainedPr[100]; /* prior on branch lengths if unconstrained */ |
|---|
| 972 | char clockPr[100]; /* prior on branch if clock enforced */ |
|---|
| 973 | char clockVarPr[100]; /* prior on clock rate variation (strict, cpp, mb(rateautocorr)) */ |
|---|
| 974 | char nodeAgePr[100]; /* prior on node depths (unconstrained, constraints) */ |
|---|
| 975 | char speciationPr[100]; /* prior on speciation rate */ |
|---|
| 976 | MrBFlt speciationFix; |
|---|
| 977 | MrBFlt speciationUni[2]; |
|---|
| 978 | MrBFlt speciationExp; |
|---|
| 979 | char extinctionPr[100]; /* prior on extinction rate */ |
|---|
| 980 | MrBFlt extinctionFix; |
|---|
| 981 | MrBFlt extinctionBeta[2]; |
|---|
| 982 | MrBFlt extinctionExp; |
|---|
| 983 | char sampleStrat[30]; /* taxon sampling strategy (for b-d process) */ |
|---|
| 984 | MrBFlt sampleProb; /* taxon sampling fraction (for b-d process) */ |
|---|
| 985 | char treeAgePr[100]; /* prior on tree age for uniform clock prior */ |
|---|
| 986 | MrBFlt treeAgeGamma[2]; |
|---|
| 987 | MrBFlt treeAgeExp; |
|---|
| 988 | MrBFlt treeAgeFix; |
|---|
| 989 | char clockRatePr[100]; /* prior on base substitution rate of tree for clock trees */ |
|---|
| 990 | MrBFlt clockRateNormal[2]; |
|---|
| 991 | MrBFlt clockRateLognormal[2]; |
|---|
| 992 | MrBFlt clockRateGamma[2]; |
|---|
| 993 | MrBFlt clockRateExp; |
|---|
| 994 | MrBFlt clockRateFix; |
|---|
| 995 | char popSizePr[100]; /* prior on population size */ |
|---|
| 996 | MrBFlt popSizeFix; |
|---|
| 997 | MrBFlt popSizeUni[2]; |
|---|
| 998 | MrBFlt popSizeLognormal[2]; |
|---|
| 999 | MrBFlt popSizeGamma[2]; |
|---|
| 1000 | MrBFlt popSizeNormal[2]; |
|---|
| 1001 | char popVarPr[100]; /* prior on pop. size variation across tree */ |
|---|
| 1002 | char growthPr[100]; /* prior on coalescence growth rate */ |
|---|
| 1003 | MrBFlt growthFix; |
|---|
| 1004 | MrBFlt growthUni[2]; |
|---|
| 1005 | MrBFlt growthExp; |
|---|
| 1006 | MrBFlt growthNorm[2]; |
|---|
| 1007 | char cppRatePr[100]; /* prior on CPP rate */ |
|---|
| 1008 | MrBFlt cppRateFix; |
|---|
| 1009 | MrBFlt cppRateExp; |
|---|
| 1010 | char cppMultDevPr[100]; /* prior on CPP rate multiplier Lognormal variance */ |
|---|
| 1011 | MrBFlt cppMultDevFix; |
|---|
| 1012 | char tk02varPr[100]; /* prior on TK02 lognormal rate variance */ |
|---|
| 1013 | MrBFlt tk02varFix; |
|---|
| 1014 | MrBFlt tk02varUni[2]; |
|---|
| 1015 | MrBFlt tk02varExp; |
|---|
| 1016 | char igrvarPr[100]; /* prior on IGR gamma distribution variance */ |
|---|
| 1017 | MrBFlt igrvarFix; |
|---|
| 1018 | MrBFlt igrvarUni[2]; |
|---|
| 1019 | MrBFlt igrvarExp; |
|---|
| 1020 | |
|---|
| 1021 | char tratioFormat[30]; /* format used to report tratio */ |
|---|
| 1022 | char revmatFormat[30]; /* format used to report revmat */ |
|---|
| 1023 | char ratemultFormat[30]; /* format used to report ratemult */ |
|---|
| 1024 | char treeFormat[30]; /* format used to report trees/topologies */ |
|---|
| 1025 | char inferAncStates[5]; /* should ancestral states be inferred (Yes/No)?*/ |
|---|
| 1026 | char inferSiteOmegas[5]; /* should site omega vals be inferred (Yes/No)? */ |
|---|
| 1027 | char inferSiteRates[5]; /* should site rates be inferred (Yes/No)? */ |
|---|
| 1028 | char inferPosSel[5]; /* should site selection be inferred (Yes/No)? */ |
|---|
| 1029 | } Model, ModelParams; |
|---|
| 1030 | |
|---|
| 1031 | typedef struct chain |
|---|
| 1032 | { |
|---|
| 1033 | int numGen; /* number of MCMC cycles */ |
|---|
| 1034 | int sampleFreq; /* frequency to sample chain */ |
|---|
| 1035 | int printFreq; /* frequency to print chain */ |
|---|
| 1036 | int swapFreq; /* frequency to attempt swap of states */ |
|---|
| 1037 | int numRuns; /* number of runs */ |
|---|
| 1038 | int numChains; /* number of chains */ |
|---|
| 1039 | int isSS; /* do we do Steppingstone Sampling */ |
|---|
| 1040 | int startFromPriorSS; /* If Yes SS is moving from Prior to Posterior */ |
|---|
| 1041 | int numStepsSS; /* Number of steps in SS */ |
|---|
| 1042 | int burninSS; /* Fixed burnin for SS */ |
|---|
| 1043 | MrBFlt alphaSS; /* Beta values are distributed according to quantiles of Beta(alphaSS,1.0) distribution */ |
|---|
| 1044 | int backupCheckSS; /* Frequency of checkpoints backup */ |
|---|
| 1045 | MrBFlt chainTemp; /* chain temperature */ |
|---|
| 1046 | int userDefinedTemps; /* should we use the users temperatures? */ |
|---|
| 1047 | MrBFlt userTemps[MAX_CHAINS]; /* user-defined chain temperatures */ |
|---|
| 1048 | char chainFileName[100]; /* chain file name for output */ |
|---|
| 1049 | int chainBurnIn; /* chain burn in length */ |
|---|
| 1050 | int numStartPerts; /* number of perturbations to starting tree */ |
|---|
| 1051 | char startTree[100]; /* starting tree for chain (current/random) */ |
|---|
| 1052 | char startParams[100]; /* starting values for chain (current/reset) */ |
|---|
| 1053 | int saveBrlens; /* should branch lengths be saved */ |
|---|
| 1054 | MrBFlt weightScheme[3]; /* percent chars to increase/decrease in weight */ |
|---|
| 1055 | int calcPbf; /* should we calculate the pseudo Bayes factor */ |
|---|
| 1056 | int pbfInitBurnin; /* initial burnin when calculating pseudo BF */ |
|---|
| 1057 | int pbfSampleFreq; /* sample frequency for pseudo BF */ |
|---|
| 1058 | int pbfSampleTime; /* how many cycles to calcualate site prob. */ |
|---|
| 1059 | int pbfSampleBurnin; /* burnin period for each site for pseudo BF */ |
|---|
| 1060 | int swapAdjacentOnly; /* whether we only swap adjacent temperatures */ |
|---|
| 1061 | int redirect; /* should output go to stdout */ |
|---|
| 1062 | int allChains; /* should stats be output for all chains? */ |
|---|
| 1063 | int numSwaps; /* number of swaps to try each time */ |
|---|
| 1064 | int mcmcDiagn; /* should mcmc diagnostics be output? */ |
|---|
| 1065 | int diagnFreq; /* mcmc diagnostics frequency */ |
|---|
| 1066 | int diagnStat; /* statistic to use for mcmc diagnostics */ |
|---|
| 1067 | int relativeBurnin; /* should a relative burnin be used ? */ |
|---|
| 1068 | MrBFlt burninFraction; /* the sample fraction to discard as burnin */ |
|---|
| 1069 | int allComps; /* top conv diagnosis for all pairs? */ |
|---|
| 1070 | MrBFlt minPartFreq; /* minimum partition frequency for conv diagn */ |
|---|
| 1071 | MrBFlt stopVal; /* top conv diagn value to reach before stopping */ |
|---|
| 1072 | int stopRule; /* use stop rule? */ |
|---|
| 1073 | STATS *stat; /* ptr to structs with mcmc diagnostics info */ |
|---|
| 1074 | Tree *dtree; /* pointing to tree used for conv diagnostics */ |
|---|
| 1075 | TreeList *treeList; /* vector of tree lists for saving trees */ |
|---|
| 1076 | int saveTrees; /* save tree samples for later removal? */ |
|---|
| 1077 | int stopTreeGen; /* generation after which no trees need be saved */ |
|---|
| 1078 | fpos_t *tFilePos; /* position for reading trees for removal */ |
|---|
| 1079 | int printMax; /* maximum number of chains to print */ |
|---|
| 1080 | int printAll; /* whether to print all or only cold chains */ |
|---|
| 1081 | int checkPoint; /* should we use check-pointing? */ |
|---|
| 1082 | int checkFreq; /* check-pointing frequency */ |
|---|
| 1083 | int runWithData; /* should we run with or without data? */ |
|---|
| 1084 | int orderTaxa; /* order taxa before printing tree to file? */ |
|---|
| 1085 | int append; /* order taxa before printing tree to file? */ |
|---|
| 1086 | int autotune; /* autotune tuning parameters of proposals ? */ |
|---|
| 1087 | int tuneFreq; /* autotuning frequency */ |
|---|
| 1088 | } Chain; |
|---|
| 1089 | |
|---|
| 1090 | typedef struct modelinfo |
|---|
| 1091 | { |
|---|
| 1092 | /* General model information */ |
|---|
| 1093 | int dataType; /* data type for partition */ |
|---|
| 1094 | int nucModelId; /* structure of nucleotide model */ |
|---|
| 1095 | int nst; /* # substitution types */ |
|---|
| 1096 | int aaModelId; /* amino acid model type */ |
|---|
| 1097 | int parsModelId; /* is parsimony model used YES/NO */ |
|---|
| 1098 | |
|---|
| 1099 | /* Specific model information */ |
|---|
| 1100 | int numGammaCats; /* number of gamma cats (1 if inapplic.) */ |
|---|
| 1101 | int numBetaCats; /* number of beta cats (1 if inapplic.) */ |
|---|
| 1102 | int numOmegaCats; /* number of omega cats (1 if inapplic.) */ |
|---|
| 1103 | int numTiCats; /* number of cats needing different tis */ |
|---|
| 1104 | int numModelStates; /* number of states including hidden ones */ |
|---|
| 1105 | int numStates; /* number of observable discrete states */ |
|---|
| 1106 | |
|---|
| 1107 | /* Model parameter pointers */ |
|---|
| 1108 | Param *tRatio; /* ptr to tRatio used in model */ |
|---|
| 1109 | Param *revMat; /* ptr to revMat used in model */ |
|---|
| 1110 | Param *omega; /* ptr to omega used in model */ |
|---|
| 1111 | Param *stateFreq; /* ptr to statFreq used in model */ |
|---|
| 1112 | Param *shape; /* ptr to shape used in model */ |
|---|
| 1113 | Param *pInvar; /* ptr to pInvar used in model */ |
|---|
| 1114 | Param *correlation; /* ptr to correlation used in model */ |
|---|
| 1115 | Param *switchRates; /* ptr to switchRates (off<->on) */ |
|---|
| 1116 | Param *rateMult; /* ptr to parition rateMult used in model */ |
|---|
| 1117 | Param *geneTreeRateMult; /* ptr to gene tree rateMult used in model */ |
|---|
| 1118 | Param *speciationRates; /* ptr to speciationRates used in model */ |
|---|
| 1119 | Param *extinctionRates; /* ptr to extinctionRates used in model */ |
|---|
| 1120 | Param *popSize; /* ptr to population size used in model */ |
|---|
| 1121 | Param *growthRate; /* ptr to growth rate used in model */ |
|---|
| 1122 | Param *topology; /* ptr to topology used in model */ |
|---|
| 1123 | Param *brlens; /* ptr to brlens (and tree) used in model */ |
|---|
| 1124 | Param *speciesTree; /* ptr to species tree used in model */ |
|---|
| 1125 | Param *aaModel; /* ptr to amino acid matrix used */ |
|---|
| 1126 | Param *cppMultDev; /* ptr to cpp ratemult lognormal variance */ |
|---|
| 1127 | Param *cppRate; /* ptr to CPP rate used in model */ |
|---|
| 1128 | Param *cppEvents; /* ptr to CPP events */ |
|---|
| 1129 | Param *tk02var; /* ptr to variance for TK02 relaxed clock */ |
|---|
| 1130 | Param *tk02BranchRates; /* ptr to branch rates for TK02 relaxed clock */ |
|---|
| 1131 | Param *igrvar; /* ptr to gamma var for IGR relaxed clock */ |
|---|
| 1132 | Param *igrBranchRates; /* ptr to branch rates for IGR relaxed clock*/ |
|---|
| 1133 | Param *clockRate; /* ptr to clock rate parameter */ |
|---|
| 1134 | |
|---|
| 1135 | /* Information about characters and transformations */ |
|---|
| 1136 | int numChars; /* number of compressed characters */ |
|---|
| 1137 | int numUncompressedChars; /* number of uncompressed characters */ |
|---|
| 1138 | int numDummyChars; /* number of dummy characters */ |
|---|
| 1139 | int compMatrixStart; /* start column in compressed matrix */ |
|---|
| 1140 | int compMatrixStop; /* stop column in compressed matrix */ |
|---|
| 1141 | int compCharStart; /* start char among compressed chars */ |
|---|
| 1142 | int compCharStop; /* stop char among compressed chars */ |
|---|
| 1143 | int parsMatrixStart; /* start column in parsimony matrix */ |
|---|
| 1144 | int parsMatrixStop; /* stop collumn in parsimony matrix */ |
|---|
| 1145 | int nParsIntsPerSite; /* # parsimony ints per character */ |
|---|
| 1146 | int nCharsPerSite; /* number chars per site (eg 3 for codon) */ |
|---|
| 1147 | int rateProbStart; /* start of rate probs (for adgamma) */ |
|---|
| 1148 | |
|---|
| 1149 | /* Variables for eigen decompositions */ |
|---|
| 1150 | int cijkLength; /* stores length of cijk vector */ |
|---|
| 1151 | int nCijkParts; /* stores number of cijk partitions (1 except for omega/covarion models) */ |
|---|
| 1152 | int upDateCijk; /* whether cijk vector needs to be updated */ |
|---|
| 1153 | |
|---|
| 1154 | /* Variables for standard model */ |
|---|
| 1155 | int *tiIndex; /* index to trans probs for each compressed char*/ |
|---|
| 1156 | int *bsIndex; /* index to stat freqs for each compressed char */ |
|---|
| 1157 | int *nStates; /* # states of each compressed char */ |
|---|
| 1158 | int *cType; /* whether char is ord, unord or irrev */ |
|---|
| 1159 | int *weight; /* prior weight of each compressed char */ |
|---|
| 1160 | int isTiNeeded[20]; /* marks whether a trans prob matrix is needed */ |
|---|
| 1161 | |
|---|
| 1162 | /* Gibbs sampling of gamma site rate parameters */ |
|---|
| 1163 | CLFlt ***catLike; /* likelihood for Gibbs sampling of gamma */ |
|---|
| 1164 | CLFlt ***catLnScaler; /* scaler for Gibbs sampling of gamma */ |
|---|
| 1165 | int gibbsGamma; /* flags whether Gibbs sampling of discrete gamma is used */ |
|---|
| 1166 | int gibbsFreq; /* frequency of Gibbs resampling of discrete gamma */ |
|---|
| 1167 | |
|---|
| 1168 | /* Variables for parsimony sets and parsimony calculations */ |
|---|
| 1169 | MrBFlt parsTreeLength[MAX_CHAINS*2]; /* parsimony tree lengths for chains */ |
|---|
| 1170 | SafeLong **parsSets; /* parsimony sets */ |
|---|
| 1171 | int numParsSets; /* number of parsimony sets */ |
|---|
| 1172 | CLFlt *parsNodeLens; /* parsimony node lengths */ |
|---|
| 1173 | int numParsNodeLens; /* number of parsimony node lengths */ |
|---|
| 1174 | |
|---|
| 1175 | /* Miscellaneous parameters */ |
|---|
| 1176 | int mark; /* scratch parameter */ |
|---|
| 1177 | int parsimonyBasedMove; /* is parsimony-based move used (YES/NO) */ |
|---|
| 1178 | |
|---|
| 1179 | /* Variables for conditional likelihoods */ |
|---|
| 1180 | int upDateCl; /* flags whether update of cond likes needed */ |
|---|
| 1181 | int upDateAll; /* flags whether update of entire tree is needed*/ |
|---|
| 1182 | int *isPartAmbig; /* is tip partially ambiguous? */ |
|---|
| 1183 | int **termState; /* index arrays for terminal branch shortcuts */ |
|---|
| 1184 | CLFlt *invCondLikes; /* space for the invariable cond likes */ |
|---|
| 1185 | CLFlt **condLikes; /* space for the cond likes */ |
|---|
| 1186 | CLFlt **tiProbs; /* space for the ti probs */ |
|---|
| 1187 | CLFlt **scalers; /* space for the node and site scalers */ |
|---|
| 1188 | CLFlt **clP; /* handy pointers to cond likes for ti cats */ |
|---|
| 1189 | #if defined (SSE_ENABLED) |
|---|
| 1190 | __m128 **clP_SSE; /* handy pointers to cond likes, SSE version */ |
|---|
| 1191 | int numSSEChars; /* number of compact SSE character groups */ |
|---|
| 1192 | CLFlt *lnL_SSE; /* temp storage for log site likes */ |
|---|
| 1193 | CLFlt *lnLI_SSE; /* temp storage for log site invariable likes */ |
|---|
| 1194 | #endif |
|---|
| 1195 | MrBFlt **cijks; /* space for cijks */ |
|---|
| 1196 | int **condLikeIndex; /* index to cond like space for nodes & chains */ |
|---|
| 1197 | int *condLikeScratchIndex; /* index to scratch space for node cond likes */ |
|---|
| 1198 | int **nodeScalerIndex; /* index to scaler space for nodes & chains */ |
|---|
| 1199 | int *nodeScalerScratchIndex; /* index to scratch space for node scalers */ |
|---|
| 1200 | int **scalersSet; /* flags whether scalers are set */ |
|---|
| 1201 | int *scalersSetScratch; /* scratch flag for whether scalers are set */ |
|---|
| 1202 | int *siteScalerIndex; /* index to site scaler space for chains */ |
|---|
| 1203 | int siteScalerScratchIndex; /* index to scratch space for site scalers */ |
|---|
| 1204 | int **tiProbsIndex; /* index to ti probs for branches & chains */ |
|---|
| 1205 | int *tiProbsScratchIndex; /* index to scratch space for branch ti probs */ |
|---|
| 1206 | int *cijkIndex; /* index to cijks for chains */ |
|---|
| 1207 | int cijkScratchIndex; /* index to scratch space for cijks */ |
|---|
| 1208 | int numCondLikes; /* number of cond like arrays */ |
|---|
| 1209 | int numScalers; /* number of scaler arrays */ |
|---|
| 1210 | int numTiProbs; /* number of ti prob arrays */ |
|---|
| 1211 | int condLikeLength; /* length of cond like array (incl. ti cats) */ |
|---|
| 1212 | int tiProbLength; /* length of ti prob array */ |
|---|
| 1213 | MrBFlt lnLike[MAX_CHAINS]; /* log like for chain */ |
|---|
| 1214 | CLFlt *ancStateCondLikes; /* ancestral state cond like array */ |
|---|
| 1215 | |
|---|
| 1216 | /* Likelihood function pointers */ |
|---|
| 1217 | LikeDownFxn CondLikeDown; /* function for calculating partials */ |
|---|
| 1218 | LikeRootFxn CondLikeRoot; /* function for calculating partials at root */ |
|---|
| 1219 | LikeScalerFxn CondLikeScaler; /* function for scaling partials */ |
|---|
| 1220 | LikeFxn Likelihood; /* function for getting cond likes for tree */ |
|---|
| 1221 | TiProbFxn TiProbs; /* function for calculating transition probs */ |
|---|
| 1222 | LikeUpFxn CondLikeUp; /* final-pass calculation of cond likes */ |
|---|
| 1223 | PrintAncStFxn PrintAncStates; /* function for sampling ancestral states */ |
|---|
| 1224 | StateCodeFxn StateCode; /* function for getting states from codes */ |
|---|
| 1225 | PrintSiteRateFxn PrintSiteRates; /* function for samling site rates */ |
|---|
| 1226 | |
|---|
| 1227 | /* Report variables */ |
|---|
| 1228 | int printAncStates; /* should ancestral states be printed (YES/NO) */ |
|---|
| 1229 | int printSiteRates; /* should site rates be printed (YES/NO) */ |
|---|
| 1230 | int printPosSel; /* should selection be printed (YES/NO) */ |
|---|
| 1231 | int printSiteOmegas; /* should site omegas be printed (YES/NO) */ |
|---|
| 1232 | |
|---|
| 1233 | /* likelihood calculator flags */ |
|---|
| 1234 | int useBeagle; /* use Beagle for this partition? */ |
|---|
| 1235 | int useSSE; /* use SSE for this partition? */ |
|---|
| 1236 | |
|---|
| 1237 | #if defined (BEAGLE_ENABLED) |
|---|
| 1238 | /* Beagle variables */ |
|---|
| 1239 | int useBeagleResource; /* try to use this BEAGLE resource number */ |
|---|
| 1240 | MrBFlt* branchLengths; /* array of branch lengths for Beagle */ |
|---|
| 1241 | MrBFlt* inRates; /* array of category rates for Beagle */ |
|---|
| 1242 | int* tiProbIndices; /* array of trans prob indices for Beagle */ |
|---|
| 1243 | MrBFlt* logLikelihoods; /* array of log likelihoods from Beagle */ |
|---|
| 1244 | int beagleInstance; /* beagle instance for division */ |
|---|
| 1245 | MrBFlt* inWeights; /* array of weights for Beagle root likelihood */ |
|---|
| 1246 | int* bufferIndices; /* array of partial indices for root likelihood */ |
|---|
| 1247 | int* eigenIndices; /* array of eigen indices for root likelihood */ |
|---|
| 1248 | int* childBufferIndices; /* array of child partial indices (unrooted) */ |
|---|
| 1249 | int* childTiProbIndices; /* array of child ti prob indices (unrooted) */ |
|---|
| 1250 | int* cumulativeScaleIndices; /* array of cumulative scale indices */ |
|---|
| 1251 | int rescaleBeagleAll; /* set to rescale all nodes */ |
|---|
| 1252 | int* rescaleFreq; /* rescale frequency for each chain's tree */ |
|---|
| 1253 | int rescaleFreqOld; /* holds rescale frequency of current state */ |
|---|
| 1254 | int recalculateScalers; /* shoud we recalculate scalers for current state YES/NO */ |
|---|
| 1255 | int* succesCount; /* count number of succesful computation since last reset of scalers */ |
|---|
| 1256 | int** isScalerNode; /* for each node and chain set to YES if scaled node */ |
|---|
| 1257 | int* isScalerNodeScratch; /* scratch space to hold isScalerNode of proposed state*/ |
|---|
| 1258 | long* beagleComputeCount; /* count of number of calls to likelihood */ |
|---|
| 1259 | #endif |
|---|
| 1260 | |
|---|
| 1261 | } ModelInfo; |
|---|
| 1262 | |
|---|
| 1263 | /* TODO: Delete these old pointers to cond likes, ti probs and scalers */ |
|---|
| 1264 | |
|---|
| 1265 | typedef struct sumt |
|---|
| 1266 | { |
|---|
| 1267 | int *absentTaxa; /* information on absent taxa */ |
|---|
| 1268 | int brlensDef; /* branch lengths defined? */ |
|---|
| 1269 | char sumtFileName[100]; /* name of input file */ |
|---|
| 1270 | char sumtOutfile[120]; /* name of output file */ |
|---|
| 1271 | char curFileName[120]; /* name of file being processed */ |
|---|
| 1272 | int burnin; /* actual burnin when parsing tree files */ |
|---|
| 1273 | char sumtConType[100]; /* consensus tree type */ |
|---|
| 1274 | int calcTreeprobs; /* should the individual tree probs be calculated*/ |
|---|
| 1275 | int showSumtTrees; /* should the individual tree probs be shown */ |
|---|
| 1276 | int printBrlensToFile; /* should branch lengths be printed to file */ |
|---|
| 1277 | MrBFlt brlensFreqDisplay; /* threshold for printing branch lengths to file */ |
|---|
| 1278 | int numRuns; /* number of independent analyses to summarize */ |
|---|
| 1279 | int numTrees; /* number of tree params to summarize */ |
|---|
| 1280 | int orderTaxa; /* order taxa in trees? */ |
|---|
| 1281 | MrBFlt minPartFreq; /* minimum part. freq. for overall diagnostics */ |
|---|
| 1282 | int table; /* show table of partition frequencies? */ |
|---|
| 1283 | int summary; /* show summary diagnostics ? */ |
|---|
| 1284 | int showConsensus; /* show consensus trees ? */ |
|---|
| 1285 | int consensusFormat; /* format of consensus tree */ |
|---|
| 1286 | PolyTree *tree; /* for storing tree read from file */ |
|---|
| 1287 | int *order; /* for storing topology read from file */ |
|---|
| 1288 | int orderLen; /* length of order array */ |
|---|
| 1289 | int numTreesInLastBlock; /* number of trees in last block */ |
|---|
| 1290 | int numTreesEncountered; /* number of trees encounted in total */ |
|---|
| 1291 | int numTreesSampled; /* number of sampled trees in total */ |
|---|
| 1292 | int isRooted; /* is sumt tree rooted ? */ |
|---|
| 1293 | int isRelaxed; /* is sumt tree a relaxed clock tree ? */ |
|---|
| 1294 | int isClock; /* is sumt tree a clock tree ? */ |
|---|
| 1295 | int isCalibrated; /* is sumt tree calibrated ? */ |
|---|
| 1296 | int nESets; /* number of event sets */ |
|---|
| 1297 | int nBSets; /* number of branch rate sets */ |
|---|
| 1298 | char **bSetName; /* name of effective branch length sets */ |
|---|
| 1299 | char **eSetName; /* name of event sets */ |
|---|
| 1300 | int popSizeSet; /* do sumt trees have population size set? */ |
|---|
| 1301 | char *popSizeSetName; /* name of population size set */ |
|---|
| 1302 | int SafeLongsNeeded; /* number of safe longs needed for taxon bits */ |
|---|
| 1303 | int runId; /* id of run being processed */ |
|---|
| 1304 | int numTaxa; /* number of sumt taxa */ |
|---|
| 1305 | char **taxaNames; /* names of sumt taxa */ |
|---|
| 1306 | int *numFileTrees; /* number of trees per file */ |
|---|
| 1307 | int *numFileTreesSampled; /* number of trees sampled per file */ |
|---|
| 1308 | int HPD; /* use highest posterior density? */ |
|---|
| 1309 | } Sumt; |
|---|
| 1310 | |
|---|
| 1311 | /* formats for consensus trees */ |
|---|
| 1312 | #define SIMPLE 0 |
|---|
| 1313 | #define FIGTREE 1 |
|---|
| 1314 | |
|---|
| 1315 | typedef struct comptree |
|---|
| 1316 | { |
|---|
| 1317 | char comptFileName1[120]; /* name of first input file */ |
|---|
| 1318 | char comptFileName2[120]; /* name of second input file */ |
|---|
| 1319 | char comptOutfile[120]; /* name of output file */ |
|---|
| 1320 | int burnin; /* actual burnin used when parsing tree files */ |
|---|
| 1321 | MrBFlt minPartFreq; /* use partitions with frequency >= minPartFreq */ |
|---|
| 1322 | } Comptree; |
|---|
| 1323 | |
|---|
| 1324 | typedef struct sump |
|---|
| 1325 | { |
|---|
| 1326 | char sumpFileName[100]; /* name of input file */ |
|---|
| 1327 | char sumpOutfile[120]; /* name of output file */ |
|---|
| 1328 | //int plot; /* output plot (y/n)? */ |
|---|
| 1329 | int table; /* output table (y/n)? */ |
|---|
| 1330 | int margLike; /* output marginal likelihood (y/n)? */ |
|---|
| 1331 | int numRuns; /* number of independent analyses to summarize */ |
|---|
| 1332 | int allRuns; /* should data for all runs be printed (yes/no)? */ |
|---|
| 1333 | int HPD; /* use highest posterior density? */ |
|---|
| 1334 | MrBFlt minProb; /* cut-off for model probabilities to show */ |
|---|
| 1335 | } Sump; |
|---|
| 1336 | |
|---|
| 1337 | typedef struct sumss |
|---|
| 1338 | { |
|---|
| 1339 | //int plot; /* output plot (y/n)? */ |
|---|
| 1340 | int numRuns; /* number of independent analyses to summarize */ |
|---|
| 1341 | int allRuns; /* should data for all runs be printed (yes/no)? */ |
|---|
| 1342 | int stepToPlot; /* Which step to plot in the step plot */ |
|---|
| 1343 | int askForMorePlots; /* Should user be asked to plot for different discardfraction (y/n)? */ |
|---|
| 1344 | int smoothing; /* An integer indicating number of neighbors to average over when dooing smoothing of curvs on plots */ |
|---|
| 1345 | MrBFlt discardFraction; /* Proportion of samples discarded when ploting step plot.*/ |
|---|
| 1346 | } Sumss; |
|---|
| 1347 | |
|---|
| 1348 | |
|---|
| 1349 | typedef struct plot |
|---|
| 1350 | { |
|---|
| 1351 | char plotFileName[120]; /* name of input file */ |
|---|
| 1352 | char parameter[100]; /* parameter(s) to be plotted */ |
|---|
| 1353 | char match[100]; /* whether the match needs to be perfect */ |
|---|
| 1354 | } Plot; |
|---|
| 1355 | |
|---|
| 1356 | typedef struct |
|---|
| 1357 | { |
|---|
| 1358 | int numTrees; /* number of trees to reassemble */ |
|---|
| 1359 | int numRuns; /* number of runs to reassemble */ |
|---|
| 1360 | } ReassembleInfo; |
|---|
| 1361 | |
|---|
| 1362 | typedef struct doublet |
|---|
| 1363 | { |
|---|
| 1364 | SafeLong first, second; |
|---|
| 1365 | } Doublet; |
|---|
| 1366 | |
|---|
| 1367 | typedef struct matrix |
|---|
| 1368 | { |
|---|
| 1369 | SafeLong *origin; |
|---|
| 1370 | int rowSize; |
|---|
| 1371 | int nRows; |
|---|
| 1372 | int column; |
|---|
| 1373 | int row; |
|---|
| 1374 | } Matrix; |
|---|
| 1375 | |
|---|
| 1376 | typedef struct charinfo |
|---|
| 1377 | { |
|---|
| 1378 | int dType; |
|---|
| 1379 | int cType; |
|---|
| 1380 | int nStates; |
|---|
| 1381 | int constant[10]; |
|---|
| 1382 | int variable; |
|---|
| 1383 | int informative; |
|---|
| 1384 | } CharInfo; |
|---|
| 1385 | |
|---|
| 1386 | typedef struct |
|---|
| 1387 | { |
|---|
| 1388 | int isExcluded; /* is the character excluded */ |
|---|
| 1389 | int numStates; /* number of observed states for the character */ |
|---|
| 1390 | int charType; /* type of character */ |
|---|
| 1391 | int isMissAmbig; /* is the character missing or ambiguous */ |
|---|
| 1392 | int ctype; /* ordering of character */ |
|---|
| 1393 | int charId; /* char ID index for doublet and codon models */ |
|---|
| 1394 | int pairsId; /* char ID for doublets */ |
|---|
| 1395 | int bigBreakAfter; /* is there a large break after this character */ |
|---|
| 1396 | } |
|---|
| 1397 | CharInformation; |
|---|
| 1398 | |
|---|
| 1399 | typedef struct |
|---|
| 1400 | { |
|---|
| 1401 | int isDeleted; /* is the taxon deleted */ |
|---|
| 1402 | int charCount; /* count holder */ |
|---|
| 1403 | } |
|---|
| 1404 | TaxaInformation; |
|---|
| 1405 | |
|---|
| 1406 | typedef struct |
|---|
| 1407 | { |
|---|
| 1408 | MrBFlt curScore; |
|---|
| 1409 | MrBFlt minScore; |
|---|
| 1410 | MrBFlt totalScore; |
|---|
| 1411 | MrBFlt stopScore; |
|---|
| 1412 | MrBFlt warp; |
|---|
| 1413 | TreeNode **leaf; |
|---|
| 1414 | TreeNode **vertex; |
|---|
| 1415 | } |
|---|
| 1416 | TreeInfo; |
|---|
| 1417 | |
|---|
| 1418 | typedef struct |
|---|
| 1419 | { |
|---|
| 1420 | int allavailable; |
|---|
| 1421 | } |
|---|
| 1422 | ShowmovesParams; |
|---|
| 1423 | |
|---|
| 1424 | typedef struct |
|---|
| 1425 | { |
|---|
| 1426 | int numNames; |
|---|
| 1427 | char **names; |
|---|
| 1428 | } |
|---|
| 1429 | NameSet; |
|---|
| 1430 | |
|---|
| 1431 | #endif |
|---|