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