src/bayes.c

   1 /*
   2  *  MrBayes 3
   3  *
   4  *  (c) 2002-2013
   5  *
   6  *  John P. Huelsenbeck
   7  *  Dept. Integrative Biology
   8  *  University of California, Berkeley
   9  *  Berkeley, CA 94720-3140
  10  *  johnh@berkeley.edu
  11  *
  12  *  Fredrik Ronquist
  13  *  Swedish Museum of Natural History
  14  *  Box 50007
  15  *  SE-10405 Stockholm, SWEDEN
  16  *  fredrik.ronquist@nrm.se
  17  *
  18  *  With important contributions by
  19  *
  20  *  Paul van der Mark (paulvdm@sc.fsu.edu)
  21  *  Maxim Teslenko (maxim.teslenko@nrm.se)
  22  *
  23  *  and by many users (run 'acknowledgments' to see more info)
  24  *
  25  * This program is free software; you can redistribute it and/or
  26  * modify it under the terms of the GNU General Public License
  27  * as published by the Free Software Foundation; either version 2
  28  * of the License, or (at your option) any later version.
  29  *
  30  * This program is distributed in the hope that it will be useful,
  31  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  32  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  33  * GNU General Public License for more details (www.gnu.org).
  34  *
  35  */
  36
  37 #include "bayes.h"
  38 #include "command.h"
  39 #include "model.h"
  40 #include "sumpt.h"
  41 #include "utils.h"
  42
  43        const char* const svnRevisionBayesC = "$Rev: 1068 $";   /* Revision keyword which is expended/updated by svn on each commit/update */
  44 extern const char* const svnRevisionBestC;
  45 extern const char* const svnRevisionCommandC;
  46 extern const char* const svnRevisionLikeliC;
  47 extern const char* const svnRevisionMbbeagleC;
  48 extern const char* const svnRevisionMcmcC;
  49 extern const char* const svnRevisionModelC;
  50 extern const char* const svnRevisionProposalC;
  51 extern const char* const svnRevisionSumptC;
  52 extern const char* const svnRevisionUtilsC;
  53
  54 #ifdef HAVE_LIBREADLINE
  55 #include <readline/readline.h>
  56 #include <readline/history.h>
  57 static char **readline_completion(const char *, int, int);
  58 #endif
  59
  60 /* local prototypes */
  61 int  CommandLine (int argc, char **argv);
  62 void GetTimeSeed (void);
  63 int  InitializeMrBayes (void);
  64 void PrintHeader (void);
  65 int  ReinitializeMrBayes (void);
  66
  67 /* global variables, declared in this file */
  68 BitsLong    bitsLongWithAllBitsSet;      /* BitsLong with all bits set, for bit ops       */
  69 ModelParams defaultModel;                /* Default model; vals set in InitializeMrBayes  */
  70 int         defTaxa;                     /* flag for whether number of taxa is known      */
  71 int         defChars;                    /* flag for whether number of chars is known     */
  72 int         defMatrix;                   /* flag for whether matrix is successfull read   */
  73 int         defPartition;                /* flag for whether character partition is read  */
  74 int         defPairs;                    /* flag for whether pairs are read               */
  75 Doublet     doublet[16];                 /* holds information on states for doublets      */
  76 int         fileNameChanged;             /* has file name been changed ?                  */
  77 RandLong    globalSeed;                  /* seed that is initialized at start up          */
  78 char        **modelIndicatorParams;      /* model indicator params                        */
  79 char        ***modelElementNames;        /* names for component models                    */
  80 int         nBitsInALong;                /* number of bits in a BitsLong                  */
  81 int         nPThreads;                   /* number of pthreads to use                     */
  82 int         numUserTrees;                /* number of defined user trees                  */
  83 int         readComment;                 /* should we read comment (looking for &) ?      */
  84 int         readWord;                    /* should we read word next ?                    */
  85 RandLong    runIDSeed;                   /* seed used only for determining run ID [stamp] */
  86 RandLong    swapSeed;                    /* seed used only for determining which to swap  */
  87 int         userLevel;                   /* user level                                    */
  88 PolyTree    *userTree[MAX_NUM_USERTREES];/* array of user trees                           */
  89 char        workingDir[100];             /* working directory                             */
  90
  91 #if defined (MPI_ENABLED)
  92 int         proc_id;                     /* process ID (0, 1, ..., num_procs-1)                        */
  93 int         num_procs;                   /* number of active processors                                */
  94 MrBFlt      myStateInfo[7];              /* likelihood/prior/heat/ran/moveInfo vals of me              */
  95 MrBFlt      partnerStateInfo[7];         /* likelihood/prior/heat/ran/moveInfo vals of partner         */
  96 #endif
  97
  98 #if defined (FAST_LOG)
  99 CLFlt       scalerValue[400];
 100 CLFlt       logValue[400];
 101 #endif
 102 /* Define to use a log lookup for 4by4 nucleotide data (actually SLOWER than normal code on intel processors) */
 103
 104
 105 int main (int argc, char *argv[])
 106 {
 107     int i;
 108
 109 #   if defined (MPI_ENABLED)
 110     int     ierror;
 111 #   endif
 112
 113 #   if defined (WIN_VERSION)
 114     HANDLE scbh;
 115     BOOL ok;
 116     DWORD lastError;
 117     COORD largestWindow;
 118     CONSOLE_SCREEN_BUFFER_INFO csbi;
 119     int currBottom;
 120     char poltmp[256];
 121
 122     scbh = GetStdHandle(STD_OUTPUT_HANDLE);
 123     GetConsoleScreenBufferInfo(scbh, &csbi);
 124     currBottom      = csbi.srWindow.Bottom;
 125     largestWindow   = GetLargestConsoleWindowSize(scbh);
 126
 127     /* Allow for screen buffer 3000 lines long and 140 characters wide */
 128     csbi.dwSize.Y = 3000;
 129     csbi.dwSize.X = 140;
 130
 131     SetConsoleScreenBufferSize(scbh, csbi.dwSize);
 132     /* Allow for maximum possible screen height */
 133     csbi.srWindow.Left      = 0; /* no change relative to current value */
 134     csbi.srWindow.Top       = 0; /* no change relative to current value */
 135     csbi.srWindow.Right     = 0; /* no change relative to current value */
 136     csbi.srWindow.Bottom    = largestWindow.Y - currBottom -10; /**/
 137     ok = SetConsoleWindowInfo(scbh, FALSE, &csbi.srWindow);
 138     if (ok == FALSE)
 139         {
 140         lastError = GetLastError();
 141         GetConsoleScreenBufferInfo(scbh, &csbi);
 142         sprintf(poltmp, "\nlastError = %d", lastError);
 143         // printf (poltmp);
 144         }
 145 #   endif
 146
 147     /* calculate the size of a long - used by bit manipulation functions */
 148     nBitsInALong = sizeof(BitsLong) * 8;
 149     for (i=0; i<nBitsInALong; i++)
 150         SetBit(i, &bitsLongWithAllBitsSet);
 151
 152 #   if defined (__MWERKS__) & defined (MAC_VERSION)
 153     /* Set up interface when using the Metrowerks compiler. This
 154        should work for either Macintosh or Windows. */
 155     SIOUXSetTitle("\pMrBayes v3.2");
 156     SIOUXSettings.fontface         = 0;  /* plain=0; bold=1 */
 157     SIOUXSettings.setupmenus       = 0;
 158     SIOUXSettings.autocloseonquit  = 1;
 159     SIOUXSettings.asktosaveonclose = 0;
 160     SIOUXSettings.rows             = 60;
 161     SIOUXSettings.columns          = 90;
 162 #   endif
 163
 164 #   if defined (MPI_ENABLED)
 165     ierror = MPI_Init(&argc, &argv);
 166     if (ierror != MPI_SUCCESS)
 167         {
 168         MrBayesPrint ("%s   Problem initializing MPI\n", spacer);
 169         exit (1);
 170         }
 171     ierror = MPI_Comm_size(MPI_COMM_WORLD, &num_procs);
 172     if (ierror != MPI_SUCCESS)
 173         {
 174         MrBayesPrint ("%s   Problem getting the number of processors\n", spacer);
 175         exit (1);
 176         }
 177     ierror = MPI_Comm_rank(MPI_COMM_WORLD, &proc_id);
 178     if (ierror != MPI_SUCCESS)
 179         {
 180         MrBayesPrint ("%s   Problem getting processors rank\n", spacer);
 181         exit (1);
 182         }
 183 #   endif
 184
 185 #   ifdef HAVE_LIBREADLINE
 186     rl_attempted_completion_function = readline_completion;
 187 #   endif
 188     /* Set up parameter table. */
 189     SetUpParms ();
 190
 191     /* initialize seed using current time */
 192     GetTimeSeed ();
 193
 194     /* Initialize the variables of the program. */
 195     InitializeMrBayes ();
 196
 197     /* Print the nifty header. */
 198     PrintHeader ();
 199
 200     /* Go to the command line, process any arguments passed to the program
 201        and then wait for input. */
 202     i = CommandLine (argc, argv);
 203
 204 #   if defined (MPI_ENABLED)
 205     MPI_Finalize();
 206 #   endif
 207
 208     if (i == ERROR)
 209         return (1);
 210     else
 211         return (0);
 212 }
 213
 214
 215 int CommandLine (int argc, char **argv)
 216 {
 217     int     i, message, nProcessedArgs;
 218     char    cmdStr[CMD_STRING_LENGTH];
 219 #   ifdef HAVE_LIBREADLINE
 220 #       ifndef MPI_ENABLED
 221     char    *cmdStrP;
 222 #       endif
 223 #   endif
 224 #   if defined (MPI_ENABLED)
 225     int     ierror;
 226 #   endif
 227
 228     for (i=0;i<CMD_STRING_LENGTH;i++) cmdStr[i]='\0';
 229
 230     /* wait for user-input commands */
 231     nProcessedArgs = 1; /* first argument is program name and needs not be processed */
 232     if (nProcessedArgs < argc)
 233         {
 234         mode = NONINTERACTIVE;  /* when a command is passed into the program, the default is to exit without listening to stdin */
 235         autoClose = YES;
 236         autoOverwrite = YES;
 237         noWarn = YES;
 238         quitOnError = YES;
 239         }
 240     for (;;)
 241         {
 242         if (nProcessedArgs < argc)
 243             {
 244             /* we are here only if a command that has been passed
 245                into the program remains to be processed */
 246             if (nProcessedArgs == 1 && (strcmp(argv[1],"-i") == 0 || strcmp(argv[1],"-I") == 0))
 247                 {
 248                 mode = INTERACTIVE;
 249                 autoClose = NO;
 250                 autoOverwrite = YES;
 251                 noWarn = NO;
 252                 quitOnError = NO;
 253                 }
 254             else
 255                 sprintf (cmdStr, "Execute %s", argv[nProcessedArgs]);
 256             nProcessedArgs++;
 257             }
 258         else
 259             {
 260             /* first check if we are in noninteractive mode and quit if so */
 261             if (mode == NONINTERACTIVE)
 262                 {
 263                 MrBayesPrint ("%s   Tasks completed, exiting program because mode is noninteractive\n", spacer);
 264                 MrBayesPrint ("%s   To return control to the command line after completion of file processing, \n", spacer);
 265                 MrBayesPrint ("%s   set mode to interactive with 'mb -i <filename>' (i is for interactive)\n", spacer);
 266                 MrBayesPrint ("%s   or use 'set mode=interactive'\n\n", spacer);
 267                 return (NO_ERROR);
 268                 }
 269             /* normally, we simply wait at the prompt for a
 270                user action */
 271 #   if defined (MPI_ENABLED)
 272             if (proc_id == 0)
 273                 {
 274                 /* do not use readline because OpenMPI does not handle it */
 275                 MrBayesPrint ("MrBayes > ");
 276                 fflush (stdin);
 277                 if (fgets (cmdStr, CMD_STRING_LENGTH - 2, stdin) == NULL)
 278                     {
 279                     if (feof(stdin))
 280                         MrBayesPrint ("%s   End of File encountered on stdin; quitting\n", spacer);
 281                     else
 282                         MrBayesPrint ("%s   Could not read command from stdin; quitting\n", spacer);
 283                     strcpy (cmdStr,"quit;\n");
 284                     }
 285                 }
 286             ierror = MPI_Bcast (&cmdStr, CMD_STRING_LENGTH, MPI_CHAR, 0, MPI_COMM_WORLD);
 287             if (ierror != MPI_SUCCESS)
 288                 {
 289                 MrBayesPrint ("%s   Problem broadcasting command string\n", spacer);
 290                 }
 291 #   else
 292 #       ifdef HAVE_LIBREADLINE
 293             cmdStrP = readline("MrBayes > ");
 294             if (cmdStrP!=NULL)
 295                     {
 296                     strncpy (cmdStr,cmdStrP,CMD_STRING_LENGTH - 2);
 297                     if (*cmdStrP)
 298                         add_history (cmdStrP);
 299                     free (cmdStrP);
 300                     }
 301             else /* fall through to if (feof(stdin))..*/
 302 #       else
 303             MrBayesPrint ("MrBayes > ");
 304             fflush (stdin);
 305             if (fgets (cmdStr, CMD_STRING_LENGTH - 2, stdin) == NULL)
 306 #       endif
 307                 {
 308                 if (feof(stdin))
 309                     MrBayesPrint ("%s   End of File encountered on stdin; quitting\n", spacer);
 310                 else
 311                     MrBayesPrint ("%s   Could not read command from stdin; quitting\n", spacer);
 312                 strcpy (cmdStr,"quit;\n");
 313                 }
 314 #   endif
 315             }
 316         i = 0;
 317         while (cmdStr[i] != '\0' && cmdStr[i] != '\n')
 318             i++;
 319         cmdStr[i++] = ';';
 320         cmdStr[i] = '\0';
 321         MrBayesPrint ("\n");
 322         if (cmdStr[0] != ';')
 323             {
 324             /* check that all characters in the string are valid */
 325             if (CheckStringValidity (cmdStr) == ERROR)
 326                 {
 327                 MrBayesPrint ("   Unknown character in command string\n\n");
 328                 }
 329             else
 330                 {
 331                 expecting = Expecting(COMMAND);
 332                 message = ParseCommand (cmdStr);
 333
 334                 if (message == NO_ERROR_QUIT)
 335                     return (NO_ERROR);
 336
 337                 if (message == ERROR && quitOnError == YES)
 338                     {
 339                     MrBayesPrint ("%s   Will exit with signal 1 (error) because quitonerror is set to yes\n", spacer);
 340                     MrBayesPrint ("%s   If you want control to be returned to the command line on error,\n", spacer);
 341                     MrBayesPrint ("%s   use 'mb -i <filename>' (i is for interactive) or use 'set quitonerror=no'\n\n", spacer);
 342                     return (ERROR);
 343                     }
 344
 345 #   if defined (MPI_ENABLED)
 346                 ierror = MPI_Barrier (MPI_COMM_WORLD);
 347                 if (ierror != MPI_SUCCESS)
 348                     {
 349                     MrBayesPrint ("%s   Problem at command barrier\n", spacer);
 350                     }
 351 #   endif
 352
 353                 MrBayesPrint ("\n");
 354                 }
 355             }
 356         }
 357 }
 358
 359
 360 #ifdef HAVE_LIBREADLINE
 361 extern char *command_generator(const char *text, int state);
 362
 363 char **readline_completion (const char *text, int start, int stop)
 364 {
 365     char **matches = (char **) NULL;
 366
 367 #   ifdef COMPLETIONMATCHES
 368     if (start == 0)
 369             matches = rl_completion_matches (text, command_generator);
 370 #   endif
 371
 372     return (matches);
 373 }
 374 #endif
 375
 376
 377 unsigned FindMaxRevision (unsigned amount, ...)
 378 {
 379     const char* cur;
 380     char tmp[20];
 381     unsigned val,i,max;
 382
 383     va_list vl;
 384     va_start(vl,amount);
 385     max=0;
 386     for (i=0;i<amount;i++)
 387         {
 388         cur=va_arg(vl,const char*);
 389         sscanf(cur,"%s %d",tmp,&val);
 390         max=(max>val)?max:val;
 391         }
 392     va_end(vl);
 393     return max;
 394 }
 395
 396
 397 void GetTimeSeed (void)
 398 {
 399     time_t      curTime;
 400
 401 #   if defined (MPI_ENABLED)
 402     int         ierror;
 403
 404     if (proc_id == 0)
 405         {
 406         curTime = time(NULL);
 407         globalSeed  = (RandLong)curTime;
 408         if (globalSeed < 0)
 409             globalSeed = -globalSeed;
 410         }
 411     ierror = MPI_Bcast(&globalSeed, 1, MPI_LONG, 0, MPI_COMM_WORLD);
 412     if (ierror != MPI_SUCCESS)
 413         {
 414         MrBayesPrint ("%s   Problem broadcasting seed\n", spacer);
 415         }
 416
 417     if (proc_id == 0)
 418         {
 419         /* Note: swapSeed will often be same as globalSeed */
 420         curTime = time(NULL);
 421         swapSeed  = (RandLong)curTime;
 422         if (swapSeed < 0)
 423             swapSeed = -swapSeed;
 424         }
 425     ierror = MPI_Bcast(&swapSeed, 1, MPI_LONG, 0, MPI_COMM_WORLD);
 426     if (ierror != MPI_SUCCESS)
 427         {
 428         MrBayesPrint ("%s   Problem broadcasting swap seed\n", spacer);
 429         }
 430
 431     if (proc_id == 0)
 432         {
 433         /* Note: runIDSeed will often be same as globalSeed */
 434         curTime = time(NULL);
 435         runIDSeed  = (RandLong)curTime;
 436         if (runIDSeed < 0)
 437             runIDSeed = -runIDSeed;
 438         }
 439     ierror = MPI_Bcast(&runIDSeed, 1, MPI_LONG, 0, MPI_COMM_WORLD);
 440     if (ierror != MPI_SUCCESS)
 441         {
 442         MrBayesPrint ("%s   Problem broadcasting run ID seed\n", spacer);
 443         }
 444
 445 #   else
 446     curTime = time(NULL);
 447     globalSeed  = (RandLong)curTime;
 448     if (globalSeed < 0)
 449         globalSeed = -globalSeed;
 450
 451     /* Note: swapSeed will often be the same as globalSeed */
 452     curTime = time(NULL);
 453     swapSeed  = (RandLong)curTime;
 454     if (swapSeed < 0)
 455         swapSeed = -swapSeed;
 456
 457     /* Note: runIDSeed will often be the same as globalSeed */
 458     curTime = time(NULL);
 459     runIDSeed  = (RandLong)curTime;
 460     if (runIDSeed < 0)
 461         runIDSeed = -runIDSeed;
 462
 463 #   endif
 464 }
 465
 466
 467 int InitializeMrBayes (void)
 468 {
 469     /* this function initializes the program; only call it at the start of execution */
 470
 471     int     i, j, growthFxn[6];
 472
 473     nBitsInALong         = sizeof(BitsLong) * 8;     /* global variable: number of bits in a BitsLong */
 474     userLevel            = STANDARD_USER;            /* default user level                            */
 475
 476     readWord             = NO;                       /* should we read a word next ?                  */
 477     readComment          = NO;                       /* should we read comments? (used by tree cmd)   */
 478     fileNameChanged      = NO;                       /* file name changed ? (used by a few commands)  */
 479     echoMB               = YES;                      /* flag used by Manual to control printing       */
 480
 481 #   if defined (MPI_ENABLED)
 482     sprintf (manFileName, "commref_mb%sp.txt", VERSION_NUMBER);  /* name of command reference file    */
 483 #   else
 484     sprintf (manFileName, "commref_mb%s.txt", VERSION_NUMBER);   /* name of command reference file    */
 485 #   endif
 486
 487     for (i=0; i<NUM_ALLOCS; i++)                     /* set allocated memory to NO                    */
 488         memAllocs[i] = NO;
 489     logToFile = NO;                                  /* should screen output be logged to a file      */
 490     strcpy(logFileName, "log.out");                  /* name of the log file                          */
 491     logFileFp = NULL;                                /* file pointer to log file                      */
 492     replaceLogFile = YES;                            /* should logfile be replace/appended to         */
 493     autoClose = NO;                                  /* set default autoclose                         */
 494     autoOverwrite = YES;                             /* set default autoOverwrite                     */
 495     noWarn = NO;                                     /* set default                                   */
 496     quitOnError = NO;                                /* set default quitOnError                       */
 497     inferAncStates = NO;                             /* set default inferAncStates                    */
 498     inferSiteOmegas = NO;                            /* set default inferSiteOmegas                   */
 499     inferSiteRates = NO;                             /* set default inferSiteRates                    */
 500     inferPosSel = NO;                                /* set default inferPosSel                       */
 501     inComment = NO;                                  /* not in comment                                */
 502     numComments = 0;                                 /* no comments encountered yet                   */
 503     mode = INTERACTIVE;                              /* set default mode                              */
 504     numOpenExeFiles = 0;                             /* no execute files open yet                     */
 505     scientific = YES;                                /* print to file using scientific format?        */
 506     precision = 6;                                   /* set default precision                         */
 507     showmovesParams.allavailable = NO;               /* do not show all available moves               */
 508     strcpy(workingDir,"");                           /* working directory                             */
 509 #   if defined (BEAGLE_ENABLED)
 510 #       if defined (WIN_VERSION)
 511     tryToUseBEAGLE = NO;                             /* try to use the BEAGLE library (NO until SSE code works in Win) */
 512 #       else
 513     tryToUseBEAGLE = NO;                             /* try to use the BEAGLE library if not Win (NO untill SSE single prec. works) */
 514 #       endif
 515     beagleScalingScheme = MB_BEAGLE_SCALE_ALWAYS;    /* use BEAGLE always scaling                     */
 516     beagleFlags = BEAGLE_FLAG_PROCESSOR_CPU;         /* default to generic CPU                        */
 517     beagleResourceNumber = 99;                       /* default to auto-resource selection            */
 518     // SSE instructions do not work in Windows environment
 519     // beagleFlags |= BEAGLE_FLAG_VECTOR_SSE;        /* default to SSE code                           */
 520     beagleResource = NULL;
 521     beagleResourceCount = 0;                         /* default has no list */
 522     beagleInstanceCount = 0;                         /* no BEAGLE instances */
 523     beagleScalingFrequency = 1000;
 524 #   endif
 525 #   if defined (THREADS_ENABLED)
 526     tryToUseThreads = NO;                            /* try to use pthread with BEAGLE library        */
 527 #   endif
 528
 529     /* set the proposal information */
 530     SetUpMoveTypes ();
 531
 532     /* Set up rates for standard amino acid models, in case we need them. */
 533     if (SetAARates () == ERROR)
 534         return (ERROR);
 535
 536     /* set up doublet information */
 537     doublet[ 0].first  = 1;   doublet[ 0].second = 1;
 538     doublet[ 1].first  = 1;   doublet[ 1].second = 2;
 539     doublet[ 2].first  = 1;   doublet[ 2].second = 4;
 540     doublet[ 3].first  = 1;   doublet[ 3].second = 8;
 541     doublet[ 4].first  = 2;   doublet[ 4].second = 1;
 542     doublet[ 5].first  = 2;   doublet[ 5].second = 2;
 543     doublet[ 6].first  = 2;   doublet[ 6].second = 4;
 544     doublet[ 7].first  = 2;   doublet[ 7].second = 8;
 545     doublet[ 8].first  = 4;   doublet[ 8].second = 1;
 546     doublet[ 9].first  = 4;   doublet[ 9].second = 2;
 547     doublet[10].first  = 4;   doublet[10].second = 4;
 548     doublet[11].first  = 4;   doublet[11].second = 8;
 549     doublet[12].first  = 8;   doublet[12].second = 1;
 550     doublet[13].first  = 8;   doublet[13].second = 2;
 551     doublet[14].first  = 8;   doublet[14].second = 4;
 552     doublet[15].first  = 8;   doublet[15].second = 8;
 553
 554 #   if defined (FAST_LOG)
 555     /* set up log table */
 556     for (i=0; i<400; i++)
 557         {
 558         scalerValue[i] = (CLFlt) ldexp (1.0, 1-i);  /* offset 1 needed to deal with scaler == 1.0 */
 559         logValue[i] = (CLFlt) log (scalerValue[i]);
 560         }
 561 #   endif
 562
 563     /* user trees */
 564     for (i=0; i<MAX_NUM_USERTREES; i++)
 565         userTree[i] = NULL;
 566
 567     /* parameter values */
 568     paramValues = NULL;
 569     intValues = NULL;
 570
 571     /* Prior model settings */
 572     defaultModel.dataType = DNA;                        /* datatype                                     */
 573     defaultModel.coding = 0;                            /* ascertainment bias                           */
 574     strcpy(defaultModel.codingString, "All");           /* ascertainment bias string                    */
 575     strcpy(defaultModel.nucModel, "4by4");              /* nucleotide model                             */
 576     strcpy(defaultModel.nst, "1");                      /* number of substitution types                 */
 577     strcpy(defaultModel.aaModelPr, "Fixed");            /* amino acid model prior                       */
 578     for (i=0; i<10; i++)
 579         defaultModel.aaModelPrProbs[i] = 0.0;
 580     strcpy(defaultModel.aaModel, "Poisson");            /* amino acid model                             */
 581     strcpy(defaultModel.parsModel, "No");               /* do not use parsimony model                   */
 582     strcpy(defaultModel.geneticCode, "Universal");      /* genetic code                                 */
 583     strcpy(defaultModel.ploidy, "Diploid");             /* ploidy level                                 */
 584     strcpy(defaultModel.omegaVar, "Equal");             /* omega variation                              */
 585     strcpy(defaultModel.ratesModel, "Equal");           /* rates across sites model                     */
 586     defaultModel.numGammaCats = 4;                      /* number of categories for gamma approximation */
 587     strcpy(defaultModel.useGibbs,"No");                 /* do not use Gibbs sampling of rate cats by default */
 588     defaultModel.gibbsFreq = 100;                       /* default Gibbs sampling frequency of rate cats*/
 589     defaultModel.numBetaCats = 5;                       /* number of categories for beta approximation  */
 590     strcpy(defaultModel.covarionModel, "No");           /* use covarion model? (yes/no)                 */
 591     strcpy(defaultModel.augmentData, "No");             /* should data be augmented                     */
 592     strcpy(defaultModel.tRatioPr, "Beta");              /* prior for ti/tv rate ratio                   */
 593     defaultModel.tRatioFix = 1.0;
 594     defaultModel.tRatioDir[0] = 1.0;
 595     defaultModel.tRatioDir[1] = 1.0;
 596     strcpy(defaultModel.revMatPr, "Dirichlet");         /* prior for GTR model (nucleotides)            */
 597     for (i=0; i<6; i++)
 598         {
 599         defaultModel.revMatFix[i] = 1.0;
 600         defaultModel.revMatDir[i] = 1.0;
 601         }
 602     defaultModel.revMatSymDir = 1.0;                    /* default prior for GTR mixed model            */
 603     strcpy (defaultModel.aaRevMatPr, "Dirichlet");      /* prior for GTR model (proteins)               */
 604     for (i=0; i<190; i++)
 605         {
 606         defaultModel.aaRevMatFix[i] = 1.0;
 607         defaultModel.aaRevMatDir[i] = 1.0;
 608         }
 609     strcpy(defaultModel.omegaPr, "Dirichlet");          /* prior for omega                              */
 610     defaultModel.omegaFix = 1.0;
 611     defaultModel.omegaDir[0] = 1.0;
 612     defaultModel.omegaDir[1] = 1.0;
 613     strcpy(defaultModel.ny98omega1pr, "Beta");          /* prior for class 1 omega (Ny98 model)         */
 614     defaultModel.ny98omega1Fixed = 0.1;
 615     defaultModel.ny98omega1Beta[0] = 1.0;
 616     defaultModel.ny98omega1Beta[1] = 1.0;
 617     strcpy(defaultModel.ny98omega3pr, "Exponential");   /* prior for class 3 omega (Ny98 model)        */
 618     defaultModel.ny98omega3Fixed = 2.0;
 619     defaultModel.ny98omega3Uni[0] = 1.0;
 620     defaultModel.ny98omega3Uni[1] = 50.0;
 621     defaultModel.ny98omega3Exp = 1.0;
 622     strcpy(defaultModel.m3omegapr, "Exponential");      /* prior for all three omegas (M3 model)        */
 623     defaultModel.m3omegaFixed[0] = 0.1;
 624     defaultModel.m3omegaFixed[1] = 1.0;
 625     defaultModel.m3omegaFixed[2] = 2.0;
 626     strcpy(defaultModel.m10betapr, "Uniform");          /* prior for omega variation (M10 model)        */
 627     strcpy(defaultModel.m10gammapr, "Uniform");
 628     defaultModel.m10betaUni[0] = 0.0;
 629     defaultModel.m10betaUni[1] = 20.0;
 630     defaultModel.m10betaExp = 1.0;
 631     defaultModel.m10betaFix[0] = 1.0;
 632     defaultModel.m10betaFix[1] = 1.0;
 633     defaultModel.m10gammaUni[0] = 0.0;
 634     defaultModel.m10gammaUni[1] = 20.0;
 635     defaultModel.m10gammaExp = 1.0;
 636     defaultModel.m10gammaFix[0] = 1.0;
 637     defaultModel.m10gammaFix[1] = 1.0;
 638     defaultModel.numM10GammaCats = 4;
 639     defaultModel.numM10BetaCats = 4;
 640     strcpy(defaultModel.codonCatFreqPr, "Dirichlet");   /* prior for selection cat frequencies         */
 641     defaultModel.codonCatFreqFix[0] = 1.0/3.0;
 642     defaultModel.codonCatFreqFix[1] = 1.0/3.0;
 643     defaultModel.codonCatFreqFix[2] = 1.0/3.0;
 644     defaultModel.codonCatDir[0] = 1.0;
 645     defaultModel.codonCatDir[1] = 1.0;
 646     defaultModel.codonCatDir[2] = 1.0;
 647     strcpy(defaultModel.stateFreqPr, "Dirichlet");      /* prior for character state frequencies        */
 648     strcpy(defaultModel.stateFreqsFixType, "Equal");
 649     for (i=0; i<200; i++)
 650         {
 651         defaultModel.stateFreqsFix[i] = 0.0;
 652         defaultModel.stateFreqsDir[i] = 1.0;
 653         }
 654     defaultModel.numDirParams = 0;
 655     strcpy(defaultModel.shapePr, "Exponential");        /* prior for gamma/lnorm shape parameter        */
 656     defaultModel.shapeFix = 0.5;
 657     defaultModel.shapeUni[0] = MIN_SHAPE_PARAM;
 658     defaultModel.shapeUni[1] = MAX_SHAPE_PARAM;
 659     defaultModel.shapeExp = 1.0;
 660     strcpy(defaultModel.pInvarPr, "Uniform");           /* prior for proportion of invariable sites     */
 661     defaultModel.pInvarFix = 0.1;
 662     defaultModel.pInvarUni[0] = 0.0;
 663     defaultModel.pInvarUni[1] = 1.0;
 664     strcpy(defaultModel.adGammaCorPr, "Uniform");       /* prior for correlation param of adGamma model */
 665     defaultModel.corrFix = 0.0;
 666     defaultModel.corrUni[0] = -1.0;
 667     defaultModel.corrUni[1] = 1.0;
 668     strcpy(defaultModel.covSwitchPr, "Uniform");        /* prior for switching rates of covarion model  */
 669     defaultModel.covswitchFix[0] = 1.0;
 670     defaultModel.covswitchFix[1] = 1.0;
 671     defaultModel.covswitchUni[0] = 0.0;
 672     defaultModel.covswitchUni[1] = 100.0;
 673     defaultModel.covswitchExp = 1.0;
 674     strcpy(defaultModel.symPiPr, "Fixed");              /* prior for pi when unidentifiable states used */
 675     defaultModel.symBetaFix = -1.0;
 676     defaultModel.symBetaUni[0] = 0.0;
 677     defaultModel.symBetaUni[1] = 20.0;
 678     defaultModel.symBetaExp = 2;
 679     strcpy(defaultModel.brownCorPr, "Fixed");           /* prior on correlation of brownian model       */
 680     defaultModel.brownCorrFix = 0.0;
 681     defaultModel.brownCorrUni[0] = -1.0;
 682     defaultModel.brownCorrUni[1] = 1.0;
 683     strcpy(defaultModel.brownScalesPr, "Gammamean");    /* prior on scales of brownian model            */
 684     defaultModel.brownScalesFix = 10.0;
 685     defaultModel.brownScalesUni[0] = 0.0;
 686     defaultModel.brownScalesUni[1] = 100.0;
 687     defaultModel.brownScalesGamma[0] = 1.0;
 688     defaultModel.brownScalesGamma[1] = 10.0;
 689     defaultModel.brownScalesGammaMean = 10.0;
 690     strcpy(defaultModel.topologyPr, "Uniform");         /* prior for tree topology                      */
 691     defaultModel.topologyFix = -1;                      /* user tree index to use for fixed topology    */
 692     defaultModel.activeConstraints = NULL;              /* which constraints are active                 */
 693     strcpy(defaultModel.brlensPr, "Unconstrained");     /* prior on branch lengths                      */
 694     defaultModel.brlensFix = -1;                        /* user tree index to use for fixed brlens      */
 695     defaultModel.brlensUni[0] = BRLENS_MIN;
 696     defaultModel.brlensUni[1] = 10.0;
 697     defaultModel.brlensExp    = 10.0;
 698     defaultModel.brlens2Exp[0]= 100.0;                  /* 1st param of twoExp prior (for internal branches) */
 699     defaultModel.brlens2Exp[1]= 10.0;                   /* 2nd param of twoExp prior (for external branches) */
 700     defaultModel.brlensDir[0] = 1.0;                    /* 1st param of GammaDir prior   */
 701 //  defaultModel.brlensDir[0] = 3.0;                    /* 1st param of invGamDir prior  */
 702     defaultModel.brlensDir[1] = 0.1;                    /* 2nd param of GammaDir prior   */
 703 //  defaultModel.brlensDir[1] = 20.0;                   /* 2nd param of invGamDir prior  */
 704     defaultModel.brlensDir[2] = 1.0;                    /* 3rd param of Dirichlet priors */
 705     defaultModel.brlensDir[3] = 1.0;                    /* 4th param of Dirichlet priors */
 706
 707     strcpy(defaultModel.unconstrainedPr, "GammaDir");   /* prior on branches if unconstrained           */
 708     strcpy(defaultModel.clockPr, "Uniform");            /* prior on branch lengths if clock enforced    */
 709     defaultModel.treeAgePr.prior = standardGamma;       /* calibration prior on tree age */
 710     strcpy(defaultModel.treeAgePr.name, "Gamma(1.00,1.00)");
 711     defaultModel.treeAgePr.priorParams[0] = 1.0;
 712     defaultModel.treeAgePr.priorParams[1] = 1.0;
 713     defaultModel.treeAgePr.priorParams[2] = -1.0;
 714     defaultModel.treeAgePr.LnPriorProb = &LnPriorProbGamma_Param_Mean_Sd;
 715     defaultModel.treeAgePr.LnPriorRatio = &LnProbRatioGamma_Param_Mean_Sd;
 716     defaultModel.treeAgePr.min = 0.0;
 717     defaultModel.treeAgePr.max = POS_INFINITY;
 718     strcpy(defaultModel.clockRatePr, "Fixed");          /* prior on base subst. rate for clock trees    */
 719     defaultModel.clockRateNormal[0] = 1.0;
 720     defaultModel.clockRateNormal[1] = 1.0;
 721     defaultModel.clockRateLognormal[0] = 0.0;           /* mean 0.0 on log scale corresponds to mean rate 1.0 */
 722     defaultModel.clockRateLognormal[1] = 0.7;           /* double or half the rate in one standard deviation  */
 723     defaultModel.clockRateGamma[0] = 1.0;
 724     defaultModel.clockRateGamma[1] = 1.0;
 725     defaultModel.clockRateExp = 1.0;
 726     defaultModel.clockRateFix = 1.0;
 727     strcpy(defaultModel.speciationPr, "Exponential");   /* prior on speciation rate (net diversification) */
 728     defaultModel.speciationFix = 0.1;
 729     defaultModel.speciationUni[0] = 0.0;
 730     defaultModel.speciationUni[1] = 10.0;
 731     defaultModel.speciationExp = 10.0;
 732     strcpy(defaultModel.extinctionPr, "Beta");          /* prior on extinction rate (turnover)          */
 733     defaultModel.extinctionFix = 0.5;
 734     defaultModel.extinctionBeta[0] = 1;
 735     defaultModel.extinctionBeta[1] = 1;
 736     strcpy(defaultModel.fossilizationPr, "Beta");       /* prior on fossilization rate (sampling proportion) */
 737     defaultModel.fossilizationFix = 0.5;
 738     defaultModel.fossilizationBeta[0] = 1;
 739     defaultModel.fossilizationBeta[1] = 1;
 740     strcpy(defaultModel.sampleStrat, "Random");         /* taxon sampling strategy                      */
 741     defaultModel.sampleProb = 1.0;                      /* extant taxon sampling fraction               */
 742     defaultModel.sampleFSNum = 0;                       /* number of fossil slice sampling events       */
 743
 744     strcpy(defaultModel.popSizePr, "Gamma");            /* prior on coalescence population size         */
 745     defaultModel.popSizeFix = 100.0;                    /* N_e = 100 */
 746     defaultModel.popSizeUni[0] = 0.0;
 747     defaultModel.popSizeUni[1] = 1000.0;
 748     defaultModel.popSizeNormal[0] = 100.0;
 749     defaultModel.popSizeNormal[1] = 30.0;
 750     defaultModel.popSizeLognormal[0] = 4.6;             /* mean on log scale corresponds to N_e = 100.0 */
 751     defaultModel.popSizeLognormal[1] = 0.4;             /* factor 10 in one standard deviation          */
 752     defaultModel.popSizeGamma[0] = 1.0;                 /* exponential with mean 1/10 = 0.1 */
 753     defaultModel.popSizeGamma[1] = 10.0;
 754     strcpy(defaultModel.popVarPr, "Equal");             /* prior on pop. size variation across tree     */
 755     strcpy(defaultModel.growthPr, "Fixed");             /* prior on coalescence growth rate prior       */
 756     defaultModel.growthFix = 0.0;
 757     defaultModel.growthUni[0] = 0.0;
 758     defaultModel.growthUni[1] = 100.0;
 759     defaultModel.growthExp = 1.0;
 760     defaultModel.growthNorm[0] = 0.0;
 761     defaultModel.growthNorm[1] = 1.0;
 762     strcpy(defaultModel.nodeAgePr, "Unconstrained");    /* prior on node depths                       */
 763     strcpy(defaultModel.clockVarPr, "Strict");          /* prior on clock rate variation              */
 764     strcpy(defaultModel.cppRatePr, "Exponential") ;     /* prior on rate of CPP for relaxed clock     */
 765     defaultModel.cppRateExp = 0.1;
 766     defaultModel.cppRateFix = 1.0;
 767     strcpy(defaultModel.cppMultDevPr, "Fixed");         /* prior on standard dev. of lognormal of rate multipliers of CPP rel clock */
 768     defaultModel.cppMultDevFix = 0.4;
 769     strcpy(defaultModel.tk02varPr, "Exponential");      /* prior on nu parameter for BM rel clock     */
 770     defaultModel.tk02varExp = 1.0;
 771     defaultModel.tk02varFix = 1.0;
 772     defaultModel.tk02varUni[0] = 0.0;
 773     defaultModel.tk02varUni[1] = 5.0;
 774     strcpy(defaultModel.igrvarPr, "Exponential");       /* prior on variance increase parameter for IGR rel clock */
 775     defaultModel.igrvarExp = 10.0;
 776     defaultModel.igrvarFix = 0.1;
 777     defaultModel.igrvarUni[0] = 0.0;
 778     defaultModel.igrvarUni[1] = 0.5;
 779     strcpy(defaultModel.mixedvarPr, "Exponential");     /* prior on var parameter for mixed rel clock */
 780     defaultModel.mixedvarExp = 10.0;
 781     defaultModel.mixedvarFix = 0.1;
 782     defaultModel.mixedvarUni[0] = 0.0;
 783     defaultModel.mixedvarUni[1] = 5.0;
 784     strcpy(defaultModel.ratePr, "Fixed");               /* prior on rate for a partition              */
 785     defaultModel.ratePrDir = 1.0;
 786     strcpy(defaultModel.generatePr, "Fixed");           /* prior on rate for a gene (multispecies coalescent) */
 787     defaultModel.generatePrDir = 1.0;
 788
 789     defaultModel.nStates = 4;                           /* number of states for partition             */
 790
 791     /* Report format settings */
 792     strcpy(defaultModel.tratioFormat, "Ratio");         /* default format for tratio                  */
 793     strcpy(defaultModel.revmatFormat, "Dirichlet");     /* default format for revmat                  */
 794     strcpy(defaultModel.ratemultFormat, "Scaled");      /* default format for ratemult                */
 795     strcpy(defaultModel.treeFormat, "Brlens");          /* default format for trees                   */
 796     strcpy(defaultModel.inferAncStates, "No");          /* do not infer ancestral states              */
 797     strcpy(defaultModel.inferPosSel, "No");             /* do not infer positive selection            */
 798     strcpy(defaultModel.inferSiteOmegas, "No");         /* do not infer site omega vals               */
 799     strcpy(defaultModel.inferSiteRates, "No");          /* do not infer site rates                    */
 800
 801     /* Allocate and initialize model indicator parameter names */
 802     modelIndicatorParams = (char **) SafeCalloc (3, sizeof (char *));
 803     modelIndicatorParams[0] = "aamodel";
 804     modelIndicatorParams[1] = "gtrsubmodel";
 805     modelIndicatorParams[2] = "";
 806
 807     /* Aamodel */
 808     modelElementNames = (char ***) SafeCalloc (3, sizeof (char **));
 809     modelElementNames[0] = (char **) SafeCalloc (12, sizeof (char *));
 810     modelElementNames[0][0]  = "Poisson";
 811     modelElementNames[0][1]  = "Jones";
 812     modelElementNames[0][2]  = "Dayhoff";
 813     modelElementNames[0][3]  = "Mtrev";
 814     modelElementNames[0][4]  = "Mtmam";
 815     modelElementNames[0][5]  = "Wag";
 816     modelElementNames[0][6]  = "Rtrev";
 817     modelElementNames[0][7]  = "Cprev";
 818     modelElementNames[0][8]  = "Vt";
 819     modelElementNames[0][9]  = "Blosum";
 820     modelElementNames[0][10] = "LG";
 821     modelElementNames[0][11] = "";
 822
 823     /* Gtrsubmodel */
 824     modelElementNames[1] = (char **) SafeCalloc (204, sizeof (char *));
 825     for (i=0; i<203; i++)
 826         {
 827         modelElementNames[1][i]  = (char *) SafeCalloc (7, sizeof (char));
 828         FromIndexToGrowthFxn(i, growthFxn);
 829         for (j=0; j<6; j++)
 830             modelElementNames[1][i][j] = '1' + growthFxn[j];
 831         modelElementNames[1][i][j] = '\0';
 832         }
 833     modelElementNames[1][203]  = "";
 834
 835     /* Termination */
 836     modelElementNames[2]    = (char **) SafeCalloc (1, sizeof(char *));
 837     modelElementNames[2][0] = "";
 838
 839     /* initialize user trees */
 840     for (i=0; i<MAX_NUM_USERTREES; i++)
 841         userTree[i] = NULL;
 842     numUserTrees = 0;
 843
 844     /* Reset translate table */
 845     ResetTranslateTable();
 846
 847     /* finally reset everything dependent on a matrix being defined */
 848     return (ReinitializeMrBayes ());
 849
 850 }
 851
 852
 853 void PrintHeader (void)
 854 {
 855     char arch[4];
 856 #   ifndef RELEASE
 857     unsigned rev = FindMaxRevision (10, svnRevisionBayesC,svnRevisionBestC,svnRevisionCommandC,svnRevisionLikeliC,svnRevisionMbbeagleC,
 858                                         svnRevisionMcmcC,svnRevisionModelC,svnRevisionProposalC,svnRevisionSumptC,svnRevisionUtilsC);
 859 #   endif
 860
 861     strcpy(arch,(sizeof(void*)==4)?"x86":"x64");
 862
 863     MrBayesPrint ("\n\n");
 864 #   ifdef RELEASE
 865     MrBayesPrint ("                            MrBayes v%s %s\n\n", VERSION_NUMBER,arch);
 866 #   else
 867     MrBayesPrint ("                        MrBayes v%s(r%d) %s\n\n", VERSION_NUMBER,rev,arch);
 868 #   endif
 869     MrBayesPrint ("                      (Bayesian Analysis of Phylogeny)\n\n");
 870 #   if defined (MPI_ENABLED)
 871     MrBayesPrint ("                             (Parallel version)\n");
 872     MrBayesPrint ("                         (%d processors available)\n\n", num_procs);
 873 #   endif
 874     MrBayesPrint ("              Distributed under the GNU General Public License\n\n\n");
 875     MrBayesPrint ("               Type \"help\" or \"help <command>\" for information\n");
 876     MrBayesPrint ("                     on the commands that are available.\n\n");
 877     MrBayesPrint ("                   Type \"about\" for authorship and general\n");
 878     MrBayesPrint ("                       information about the program.\n\n\n");
 879 }
 880
 881
 882 int ReinitializeMrBayes (void)
 883 {
 884     /* this function resets everything dependent on a matrix */
 885
 886     int             i;
 887
 888     /* reinitialize indentation */
 889     strcpy (spacer, "");                             /* holds blanks for indentation                  */
 890
 891     /* reset all taxa flags */
 892     ResetTaxaFlags();
 893
 894     /* reset all characters flags */
 895     ResetCharacterFlags();
 896
 897     /* chain parameters */
 898     chainParams.numGen = 1000000;                    /* number of MCMC cycles                         */
 899     chainParams.sampleFreq = 500;                    /* frequency to sample chain                     */
 900     chainParams.printFreq = 1000;                    /* frequency to print chain                      */
 901     chainParams.swapFreq = 1;                        /* frequency of attempting swap of states        */
 902     chainParams.numSwaps = 1;                        /* number of swaps to try each time              */
 903     chainParams.isSS = NO;
 904     chainParams.startFromPriorSS = NO;
 905     chainParams.numStepsSS = 50;
 906     chainParams.burninSS = -1;
 907     chainParams.alphaSS = 0.4;
 908     chainParams.backupCheckSS = 0;
 909     chainParams.mcmcDiagn = YES;                     /* write MCMC diagnostics to file ?              */
 910     chainParams.diagnFreq = 5000;                    /* diagnostics frequency                         */
 911     chainParams.minPartFreq = 0.10;                  /* min partition frequency for diagnostics       */
 912     chainParams.allChains = NO;                      /* calculate diagnostics for all chains ?        */
 913     chainParams.allComps = NO;                       /* do not calc diagn for all run comparisons     */
 914     chainParams.relativeBurnin = YES;                /* use relative burnin?                          */
 915     chainParams.burninFraction = 0.25;               /* default burnin fraction                       */
 916     chainParams.stopRule = NO;                       /* should stopping rule be used?                 */
 917     chainParams.stopVal = 0.05;                      /* convergence diagnostic value to reach         */
 918     chainParams.numRuns = 2;                         /* number of runs                                */
 919     chainParams.numChains = 4;                       /* number of chains                              */
 920     chainParams.chainTemp = 0.1;                     /* chain temperature                             */
 921     chainParams.redirect = NO;                       /* should printf be to stdout                    */
 922     strcpy(chainParams.chainFileName, "temp");       /* chain file name for output                    */
 923     chainParams.chainBurnIn = 0;                     /* chain burn in length                          */
 924     chainParams.numStartPerts = 0;                   /* number of perturbations to starting tree      */
 925     strcpy(chainParams.startTree, "Current");        /* starting tree for chain (random/current)      */
 926     strcpy(chainParams.startParams, "Current");      /* starting params for chain (reset/current)     */
 927     chainParams.saveBrlens = YES;                    /* should branch lengths be saved                */
 928     chainParams.weightScheme[0] = 0.0;               /* percent chars to decrease in weight           */
 929     chainParams.weightScheme[1] = 0.0;               /* percent chars to increase in weight           */
 930     chainParams.weightScheme[2] = 1.0;               /* weight increment                              */
 931     chainParams.calcPbf = NO;                        /* should we calculate the pseudo-BF?            */
 932     chainParams.pbfInitBurnin = 100000;              /* initial burnin for pseudo BF                  */
 933     chainParams.pbfSampleFreq = 10;                  /* sample frequency for pseudo BF                */
 934     chainParams.pbfSampleTime = 2000;                /* how many cycles to calcualate site prob.      */
 935     chainParams.pbfSampleBurnin = 2000;              /* burnin period for each site for pseudo BF     */
 936     chainParams.userDefinedTemps = NO;               /* should we use the users temperatures?         */
 937     for (i=0; i<MAX_CHAINS; i++)
 938         chainParams.userTemps[i] = 1.0;              /* user-defined chain temperatures               */
 939     chainParams.swapAdjacentOnly = NO;               /* swap only adjacent temperatures               */
 940     chainParams.printMax = 8;                        /* maximum number of chains to print to screen   */
 941     chainParams.printAll = YES;                      /* whether to print heated chains                */
 942     chainParams.treeList = NULL;                     /* vector of tree lists for saving trees         */
 943     chainParams.saveTrees = NO;                      /* save tree samples for later removal?          */
 944     chainParams.tFilePos = NULL;                     /* position for reading trees for removal        */
 945     chainParams.runWithData = YES;                   /* whether to run with data                      */
 946     chainParams.orderTaxa = NO;                      /* should taxa be ordered in output trees?       */
 947     chainParams.append = NO;                         /* append to previous analysis?                  */
 948     chainParams.autotune = YES;                      /* autotune?                                     */
 949     chainParams.tuneFreq = 100;                      /* autotuning frequency                          */
 950     chainParams.checkPoint = YES;                    /* should we checkpoint the run?                 */
 951     chainParams.checkFreq = 2000;                    /* check-pointing frequency                      */
 952     chainParams.diagnStat = AVGSTDDEV;               /* mcmc diagnostic to use                        */
 953
 954     /* sumt parameters */
 955     strcpy(sumtParams.sumtFileName, "temp");         /* input name for sumt command                   */
 956     strcpy(sumtParams.sumtConType, "Halfcompat");    /* type of consensus tree output                 */
 957     sumtParams.calcTreeprobs = YES;                  /* should individual tree probs be calculated    */
 958     sumtParams.showSumtTrees = NO;                   /* should the individual tree probs be shown     */
 959     sumtParams.printBrlensToFile = NO;               /* should brlens be printed to file              */
 960     sumtParams.brlensFreqDisplay = 0.50;             /* threshold for printing brlens to file         */
 961     sumtParams.numTrees = 1;                         /* number of trees to summarize                  */
 962     sumtParams.numRuns = 2;                          /* number of analyses to summarize               */
 963     sumtParams.orderTaxa = YES;                      /* order taxa in trees ?                         */
 964     sumtParams.minPartFreq = 0.10;                   /* minimum part. freq. for overall diagnostics   */
 965     sumtParams.table = YES;                          /* display table of part. freq.?                 */
 966     sumtParams.summary = YES;                        /* display overall diagnostics?                  */
 967     sumtParams.showConsensus = YES;                  /* display consensus tree(s)?                    */
 968     sumtParams.consensusFormat = FIGTREE;            /* format of consensus tree                      */
 969     strcpy (sumtParams.sumtOutfile, "temp");         /* output name for sumt command                  */
 970     sumtParams.HPD = YES;                            /* use Highest Posterior Density?                */
 971
 972     /* sump parameters */
 973     strcpy(sumpParams.sumpFileName, "temp");         /* input name for sump command                   */
 974     strcpy (sumpParams.sumpOutfile, "temp");         /* output name for sump command                  */
 975     sumpParams.numRuns = 2;                          /* number of analyses to summarize               */
 976     sumpParams.HPD = YES;                            /* use Highest Posterior Density?                */
 977     sumpParams.minProb = 0.05;                       /* min. prob. of models to include in summary    */
 978
 979     /* sumss parameters */
 980     sumssParams.numRuns= 2;                          /* number of independent analyses to summarize   */
 981     sumssParams.allRuns = YES;                       /* should data for all runs be printed (yes/no)? */
 982     sumssParams.stepToPlot = 0;                      /* Which step to plot in the step plot, 0 means burnin */
 983     sumssParams.askForMorePlots = YES;               /* Should user be asked to plot for different discardfraction (y/n)?  */
 984     sumssParams.discardFraction = 0.8;               /* Proportion of samples discarded when ploting step plot.*/
 985     sumssParams.smoothing = 0;                       /* An integer indicating number of neighbors to average over
 986                                                         when dooing smoothing of curvs on plots */
 987     /* comparetree parameters */
 988     strcpy(comptreeParams.comptFileName1, "temp.t"); /* input name for comparetree command            */
 989     strcpy(comptreeParams.comptFileName2, "temp.t"); /* input name for comparetree command            */
 990     strcpy(comptreeParams.comptOutfile, "temp.comp");/* output name for comparetree command           */
 991     comptreeParams.minPartFreq = 0.0;                /* minimum frequency of partitions to include    */
 992
 993     /* plot parameters */
 994     strcpy(plotParams.plotFileName, "temp.p");       /* input name for plot command                   */
 995     strcpy(plotParams.parameter, "lnL");             /* plotted parameter plot command                */
 996     strcpy(plotParams.match, "Perfect");             /* matching for plot command                     */
 997
 998     return (NO_ERROR);
 999 }
1000