source: branches/items/GDE/MrBAYES/mrbayes_3.2.1/mb.h

Last change on this file was 10416, checked in by aboeckma, 11 years ago

Added mr bayes (no menu yet)

File size: 61.9 KB
Line 
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
41typedef int SafeLong;
42#else
43typedef long SafeLong;
44#endif
45
46
47typedef 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 */
51typedef 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
392typedef void * VoidPtr;
393typedef int (*CmdFxn)(void);
394typedef int (*ParmFxn)(char *, char *);
395
396typedef 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 */
407enum CALPRIOR
408        {
409        unconstrained,
410        fixed,
411        offsetExponential,
412        uniform
413        };
414
415enum ConstraintType
416    {
417    PARTIAL,
418    NEGATIVE,
419    HARD
420    };
421
422/* typedef for calibration */
423typedef 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 */
436typedef struct element {
437        struct element *next;
438        int                             *order;
439} TreeListElement;
440
441/* typedef for list of trees (topologies) */
442typedef struct {
443        TreeListElement *first;
444        TreeListElement *last;
445} TreeList;
446
447/* typedef for packed tree */
448typedef 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 */
455typedef 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 */
479typedef 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 */
506typedef 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 */
528typedef 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 */
560typedef MrBFlt (*LnPriorProbFxn)(MrBFlt val, MrBFlt *priorParams);
561
562/* typedef for a ln prior prob ratio fxn */
563typedef MrBFlt (*LnPriorRatioFxn)(MrBFlt newVal, MrBFlt oldVal, MrBFlt *priorParams);
564
565/* struct for holding model parameter info for the mcmc run */
566typedef 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
614typedef 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 */
772typedef int (MoveFxn)(Param *param, int chain, SafeLong *seed, MrBFlt *lnPriorRatio, MrBFlt *lnProposalRatio, MrBFlt *mvp);
773
774/* typedef for an ApplicFxn */
775typedef int (ApplicFxn)(Param *param);
776
777/* typedef for an AutotuneFxn */
778typedef 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 */
781typedef 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 */
813typedef 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
831typedef int (*LikeDownFxn)(TreeNode *, int, int);
832typedef int (*LikeRootFxn)(TreeNode *, int, int);
833typedef int (*LikeScalerFxn)(TreeNode *, int, int);
834typedef int (*LikeFxn)(TreeNode *, int, int, MrBFlt *, int);
835typedef int (*TiProbFxn)(TreeNode *, int, int);
836typedef int (*LikeUpFxn)(TreeNode *, int, int);
837typedef int (*PrintAncStFxn)(TreeNode *, int, int);
838typedef int (*StateCodeFxn) (int);
839typedef int (*PrintSiteRateFxn) (TreeNode *, int, int);
840
841
842typedef 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       
857typedef 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
865typedef 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
1031typedef 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
1090typedef 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
1265typedef 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
1315typedef 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
1324typedef 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
1337typedef 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
1349typedef 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
1356typedef struct
1357        {
1358        int                 numTrees;                  /* number of trees to reassemble                 */
1359        int                 numRuns;                   /* number of runs to reassemble                  */
1360        } ReassembleInfo;
1361
1362typedef struct doublet
1363        {
1364        SafeLong        first, second;
1365        } Doublet;
1366
1367typedef struct matrix
1368        {
1369        SafeLong *origin;
1370        int rowSize;
1371        int nRows;
1372        int column;
1373        int row;
1374        } Matrix;
1375
1376typedef 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       
1386typedef 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
1399typedef struct 
1400        {
1401        int                     isDeleted;             /* is the taxon deleted                          */
1402        int                     charCount;             /* count holder                                  */
1403        }
1404        TaxaInformation;
1405
1406typedef 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
1418typedef struct
1419        {
1420        int             allavailable;
1421        }
1422        ShowmovesParams;
1423
1424typedef struct
1425    {
1426    int     numNames;
1427    char    **names;
1428    }
1429    NameSet;
1430
1431#endif
Note: See TracBrowser for help on using the repository browser.