7 * Dept. Integrative Biology
8 * University of California, Berkeley
9 * Berkeley, CA 94720-3140
13 * Swedish Museum of Natural History
15 * SE-10405 Stockholm, SWEDEN
16 * fredrik.ronquist@nrm.se
18 * With important contributions by
20 * Paul van der Mark (paulvdm@sc.fsu.edu)
21 * Maxim Teslenko (maxim.teslenko@nrm.se)
23 * and by many users (run 'acknowledgments' to see more info)
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.
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).
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;
54 #ifdef HAVE_LIBREADLINE
55 #include <readline/readline.h>
56 #include <readline/history.h>
57 static char **readline_completion(const char *, int, int);
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);
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 */
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 */
98 #if defined (FAST_LOG)
99 CLFlt scalerValue[400];
102 /* Define to use a log lookup for 4by4 nucleotide data (actually SLOWER than normal code on intel processors) */
105 int main (int argc, char *argv[])
109 # if defined (MPI_ENABLED)
113 # if defined (WIN_VERSION)
118 CONSOLE_SCREEN_BUFFER_INFO csbi;
122 scbh = GetStdHandle(STD_OUTPUT_HANDLE);
123 GetConsoleScreenBufferInfo(scbh, &csbi);
124 currBottom = csbi.srWindow.Bottom;
125 largestWindow = GetLargestConsoleWindowSize(scbh);
127 /* Allow for screen buffer 3000 lines long and 140 characters wide */
128 csbi.dwSize.Y = 3000;
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);
140 lastError = GetLastError();
141 GetConsoleScreenBufferInfo(scbh, &csbi);
142 sprintf(poltmp, "\nlastError = %d", lastError);
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);
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;
164 # if defined (MPI_ENABLED)
165 ierror = MPI_Init(&argc, &argv);
166 if (ierror != MPI_SUCCESS)
168 MrBayesPrint ("%s Problem initializing MPI\n", spacer);
171 ierror = MPI_Comm_size(MPI_COMM_WORLD, &num_procs);
172 if (ierror != MPI_SUCCESS)
174 MrBayesPrint ("%s Problem getting the number of processors\n", spacer);
177 ierror = MPI_Comm_rank(MPI_COMM_WORLD, &proc_id);
178 if (ierror != MPI_SUCCESS)
180 MrBayesPrint ("%s Problem getting processors rank\n", spacer);
185 # ifdef HAVE_LIBREADLINE
186 rl_attempted_completion_function = readline_completion;
188 /* Set up parameter table. */
191 /* initialize seed using current time */
194 /* Initialize the variables of the program. */
195 InitializeMrBayes ();
197 /* Print the nifty header. */
200 /* Go to the command line, process any arguments passed to the program
201 and then wait for input. */
202 i = CommandLine (argc, argv);
204 # if defined (MPI_ENABLED)
215 int CommandLine (int argc, char **argv)
217 int i, message, nProcessedArgs;
218 char cmdStr[CMD_STRING_LENGTH];
219 # ifdef HAVE_LIBREADLINE
224 # if defined (MPI_ENABLED)
228 for (i=0;i<CMD_STRING_LENGTH;i++) cmdStr[i]='\0';
230 /* wait for user-input commands */
231 nProcessedArgs = 1; /* first argument is program name and needs not be processed */
232 if (nProcessedArgs < argc)
234 mode = NONINTERACTIVE; /* when a command is passed into the program, the default is to exit without listening to stdin */
242 if (nProcessedArgs < argc)
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))
255 sprintf (cmdStr, "Execute %s", argv[nProcessedArgs]);
260 /* first check if we are in noninteractive mode and quit if so */
261 if (mode == NONINTERACTIVE)
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);
269 /* normally, we simply wait at the prompt for a
271 # if defined (MPI_ENABLED)
274 /* do not use readline because OpenMPI does not handle it */
275 MrBayesPrint ("MrBayes > ");
277 if (fgets (cmdStr, CMD_STRING_LENGTH - 2, stdin) == NULL)
280 MrBayesPrint ("%s End of File encountered on stdin; quitting\n", spacer);
282 MrBayesPrint ("%s Could not read command from stdin; quitting\n", spacer);
283 strcpy (cmdStr,"quit;\n");
286 ierror = MPI_Bcast (&cmdStr, CMD_STRING_LENGTH, MPI_CHAR, 0, MPI_COMM_WORLD);
287 if (ierror != MPI_SUCCESS)
289 MrBayesPrint ("%s Problem broadcasting command string\n", spacer);
292 # ifdef HAVE_LIBREADLINE
293 cmdStrP = readline("MrBayes > ");
296 strncpy (cmdStr,cmdStrP,CMD_STRING_LENGTH - 2);
298 add_history (cmdStrP);
301 else /* fall through to if (feof(stdin))..*/
303 MrBayesPrint ("MrBayes > ");
305 if (fgets (cmdStr, CMD_STRING_LENGTH - 2, stdin) == NULL)
309 MrBayesPrint ("%s End of File encountered on stdin; quitting\n", spacer);
311 MrBayesPrint ("%s Could not read command from stdin; quitting\n", spacer);
312 strcpy (cmdStr,"quit;\n");
317 while (cmdStr[i] != '\0' && cmdStr[i] != '\n')
322 if (cmdStr[0] != ';')
324 /* check that all characters in the string are valid */
325 if (CheckStringValidity (cmdStr) == ERROR)
327 MrBayesPrint (" Unknown character in command string\n\n");
331 expecting = Expecting(COMMAND);
332 message = ParseCommand (cmdStr);
334 if (message == NO_ERROR_QUIT)
337 if (message == ERROR && quitOnError == YES)
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);
345 # if defined (MPI_ENABLED)
346 ierror = MPI_Barrier (MPI_COMM_WORLD);
347 if (ierror != MPI_SUCCESS)
349 MrBayesPrint ("%s Problem at command barrier\n", spacer);
360 #ifdef HAVE_LIBREADLINE
361 extern char *command_generator(const char *text, int state);
363 char **readline_completion (const char *text, int start, int stop)
365 char **matches = (char **) NULL;
367 # ifdef COMPLETIONMATCHES
369 matches = rl_completion_matches (text, command_generator);
377 unsigned FindMaxRevision (unsigned amount, ...)
386 for (i=0;i<amount;i++)
388 cur=va_arg(vl,const char*);
389 sscanf(cur,"%s %d",tmp,&val);
390 max=(max>val)?max:val;
397 void GetTimeSeed (void)
401 # if defined (MPI_ENABLED)
406 curTime = time(NULL);
407 globalSeed = (RandLong)curTime;
409 globalSeed = -globalSeed;
411 ierror = MPI_Bcast(&globalSeed, 1, MPI_LONG, 0, MPI_COMM_WORLD);
412 if (ierror != MPI_SUCCESS)
414 MrBayesPrint ("%s Problem broadcasting seed\n", spacer);
419 /* Note: swapSeed will often be same as globalSeed */
420 curTime = time(NULL);
421 swapSeed = (RandLong)curTime;
423 swapSeed = -swapSeed;
425 ierror = MPI_Bcast(&swapSeed, 1, MPI_LONG, 0, MPI_COMM_WORLD);
426 if (ierror != MPI_SUCCESS)
428 MrBayesPrint ("%s Problem broadcasting swap seed\n", spacer);
433 /* Note: runIDSeed will often be same as globalSeed */
434 curTime = time(NULL);
435 runIDSeed = (RandLong)curTime;
437 runIDSeed = -runIDSeed;
439 ierror = MPI_Bcast(&runIDSeed, 1, MPI_LONG, 0, MPI_COMM_WORLD);
440 if (ierror != MPI_SUCCESS)
442 MrBayesPrint ("%s Problem broadcasting run ID seed\n", spacer);
446 curTime = time(NULL);
447 globalSeed = (RandLong)curTime;
449 globalSeed = -globalSeed;
451 /* Note: swapSeed will often be the same as globalSeed */
452 curTime = time(NULL);
453 swapSeed = (RandLong)curTime;
455 swapSeed = -swapSeed;
457 /* Note: runIDSeed will often be the same as globalSeed */
458 curTime = time(NULL);
459 runIDSeed = (RandLong)curTime;
461 runIDSeed = -runIDSeed;
467 int InitializeMrBayes (void)
469 /* this function initializes the program; only call it at the start of execution */
471 int i, j, growthFxn[6];
473 nBitsInALong = sizeof(BitsLong) * 8; /* global variable: number of bits in a BitsLong */
474 userLevel = STANDARD_USER; /* default user level */
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 */
481 # if defined (MPI_ENABLED)
482 sprintf (manFileName, "commref_mb%sp.txt", VERSION_NUMBER); /* name of command reference file */
484 sprintf (manFileName, "commref_mb%s.txt", VERSION_NUMBER); /* name of command reference file */
487 for (i=0; i<NUM_ALLOCS; i++) /* set allocated memory to 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) */
513 tryToUseBEAGLE = NO; /* try to use the BEAGLE library if not Win (NO untill SSE single prec. works) */
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;
525 # if defined (THREADS_ENABLED)
526 tryToUseThreads = NO; /* try to use pthread with BEAGLE library */
529 /* set the proposal information */
532 /* Set up rates for standard amino acid models, in case we need them. */
533 if (SetAARates () == ERROR)
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;
554 # if defined (FAST_LOG)
555 /* set up log table */
556 for (i=0; i<400; i++)
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]);
564 for (i=0; i<MAX_NUM_USERTREES; i++)
567 /* parameter values */
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 */
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) */
599 defaultModel.revMatFix[i] = 1.0;
600 defaultModel.revMatDir[i] = 1.0;
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++)
606 defaultModel.aaRevMatFix[i] = 1.0;
607 defaultModel.aaRevMatDir[i] = 1.0;
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++)
651 defaultModel.stateFreqsFix[i] = 0.0;
652 defaultModel.stateFreqsDir[i] = 1.0;
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 */
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 */
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;
789 defaultModel.nStates = 4; /* number of states for partition */
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 */
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] = "";
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] = "";
824 modelElementNames[1] = (char **) SafeCalloc (204, sizeof (char *));
825 for (i=0; i<203; i++)
827 modelElementNames[1][i] = (char *) SafeCalloc (7, sizeof (char));
828 FromIndexToGrowthFxn(i, growthFxn);
830 modelElementNames[1][i][j] = '1' + growthFxn[j];
831 modelElementNames[1][i][j] = '\0';
833 modelElementNames[1][203] = "";
836 modelElementNames[2] = (char **) SafeCalloc (1, sizeof(char *));
837 modelElementNames[2][0] = "";
839 /* initialize user trees */
840 for (i=0; i<MAX_NUM_USERTREES; i++)
844 /* Reset translate table */
845 ResetTranslateTable();
847 /* finally reset everything dependent on a matrix being defined */
848 return (ReinitializeMrBayes ());
853 void PrintHeader (void)
857 unsigned rev = FindMaxRevision (10, svnRevisionBayesC,svnRevisionBestC,svnRevisionCommandC,svnRevisionLikeliC,svnRevisionMbbeagleC,
858 svnRevisionMcmcC,svnRevisionModelC,svnRevisionProposalC,svnRevisionSumptC,svnRevisionUtilsC);
861 strcpy(arch,(sizeof(void*)==4)?"x86":"x64");
863 MrBayesPrint ("\n\n");
865 MrBayesPrint (" MrBayes v%s %s\n\n", VERSION_NUMBER,arch);
867 MrBayesPrint (" MrBayes v%s(r%d) %s\n\n", VERSION_NUMBER,rev,arch);
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);
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");
882 int ReinitializeMrBayes (void)
884 /* this function resets everything dependent on a matrix */
888 /* reinitialize indentation */
889 strcpy (spacer, ""); /* holds blanks for indentation */
891 /* reset all taxa flags */
894 /* reset all characters flags */
895 ResetCharacterFlags();
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 */
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? */
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 */
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 */
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 */