Context Navigation

bayes.c

Visit:

Last change on this file was 10416, checked in by aboeckma, 11 years ago
Added mr bayes (no menu yet)
File size: 50.4 KB

Line
1	/*
2	* MrBayes 3
3	*
4	* (c) 2002-2010
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 'acknowledgements' 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 <stdio.h>
38	#include <stdlib.h>
39	#include <time.h>
40	#include <math.h>
41	#include <string.h>
42	#include "mb.h"
43	#include "globals.h"
44	#include "bayes.h"
45	#include "command.h"
46	#include "mcmc.h"
47	#include "model.h"
48	#include "sumt.h"
49	#include "utils.h"
50
51	const char* const svnRevisionBayesC="$Rev: 502 $"; /* Revision keyword which is expended/updated by svn on each commit/update*/
52	extern const char* const svnRevisionBestC;
53	extern const char* const svnRevisionCommandC;
54	extern const char* const svnRevisionMbC;
55	extern const char* const svnRevisionMbbeagleC;
56	extern const char* const svnRevisionMbmathC;
57	extern const char* const svnRevisionMcmcC;
58	extern const char* const svnRevisionModelC;
59	extern const char* const svnRevisionPlotC;
60	extern const char* const svnRevisionSumpC;
61	extern const char* const svnRevisionSumtC;
62	extern const char* const svnRevisionTreeC;
63	extern const char* const svnRevisionUtilsC;
64
65
66	#ifdef HAVE_LIBREADLINE
67	#include <readline/readline.h>
68	#include <readline/history.h>
69	static char *readline_completion(const char , int, int);
70	#endif
71
72	#if defined (WIN_VERSION)
73	# undef NO_ERROR
74	# undef ERROR
75	# include <windows.h>
76	# include <winbase.h>
77	# undef NO_ERROR
78	# undef ERROR
79	# define NO_ERROR 0
80	# define ERROR 1
81	#endif
82
83	/* local prototypes */
84	int CommandLine (int argc, char **argv);
85	void PrintHeader (void);
86
87	/* global variables, declared in this file */
88	ModelParams defaultModel; /* Default model; values set in InitializeMrBayes */
89	int defTaxa; /* flag for whether number of taxa is known */
90	int defChars; /* flag for whether number of chars is known */
91	int defMatrix; /* flag for whether matrix is successfull read */
92	int defPartition; /* flag for whether character partition is read */
93	int defPairs; /* flag for whether pairs are read */
94	Doublet doublet[16]; /* holds information on states for doublets */
95	int fileNameChanged; /* has file name been changed ? */
96	SafeLong globalSeed; /* seed that is initialized at start up */
97	char *modelIndicatorParams; / model indicator params */
98	char **modelElementNames; / names for component models */
99	int nBitsInALong; /* number of bits in a SafeLong */
100	int nPThreads; /* number of pthreads to use */
101	int numUserTrees; /* number of defined user trees */
102	int readComment; /* should we read comment (looking for &) ? */
103	int readWord; /* should we read word next ? */
104	SafeLong runIDSeed; /* seed used only for determining run ID [stamp] */
105	SafeLong safeLongWithAllBitsSet; /* SafeLong with all bits set, for bit ops */
106	SafeLong swapSeed; /* seed used only for determining which to swap */
107	int userLevel; /* user level */
108	PolyTree userTree[MAX_NUM_USERTREES];/ array of user trees */
109	char workingDir[100]; /* working directory */
110
111	# if defined (MPI_ENABLED)
112	int proc_id; /* process ID (0, 1, ..., num_procs-1) */
113	int num_procs; /* number of active processors */
114	MrBFlt myStateInfo[7]; /* likelihood/prior/heat/ran/moveInfo vals of me */
115	MrBFlt partnerStateInfo[7]; /* likelihood/prior/heat/ran/moveInfo vals of partner */
116	# endif
117
118	#if defined (FAST_LOG)
119	CLFlt scalerValue[400];
120	CLFlt logValue[400];
121	#endif
122
123
124
125	int main (int argc, char *argv[])
126
127	{
128	int i;
129
130	# if defined (MPI_ENABLED)
131	int ierror;
132	# endif
133
134	# if defined (WIN_VERSION)
135	HANDLE scbh;
136	BOOL ok;
137	DWORD lastError;
138	COORD largestWindow;
139	CONSOLE_SCREEN_BUFFER_INFO csbi;
140	int currBottom;
141	char poltmp[256];
142
143	scbh = GetStdHandle(STD_OUTPUT_HANDLE);
144	GetConsoleScreenBufferInfo(scbh, &csbi);
145	currBottom = csbi.srWindow.Bottom;
146	largestWindow = GetLargestConsoleWindowSize(scbh);
147
148	/* Allow for screen buffer 3000 lines long and 140 characters wide */
149	csbi.dwSize.Y = 3000;
150	csbi.dwSize.X = 140;
151
152	SetConsoleScreenBufferSize(scbh, csbi.dwSize);
153	/* Allow for maximum possible screen height */
154	csbi.srWindow.Left = 0; /* no change relative to current value */
155	csbi.srWindow.Top = 0; /* no change relative to current value */
156	csbi.srWindow.Right = 0; /* no change relative to current value */
157	csbi.srWindow.Bottom = largestWindow.Y - currBottom -10; /**/
158	ok = SetConsoleWindowInfo(scbh, FALSE, &csbi.srWindow);
159	if (ok == FALSE)
160	{
161	lastError = GetLastError();
162	GetConsoleScreenBufferInfo(scbh, &csbi);
163	sprintf(poltmp, "\nlastError = %d", lastError);
164	printf(poltmp);
165	}
166	# endif
167
168	/mtrace();/
169	/* calculate the size of a long - used by bit manipulation functions */
170	nBitsInALong = sizeof(SafeLong) * 8;
171	if (nBitsInALong > 32) /* Do not use more than 32 bits until we */
172	nBitsInALong = 32; /* understand how 64-bit longs are handled. */
173	for (i=0; i<nBitsInALong; i++)
174	SetBit(i, &safeLongWithAllBitsSet);
175
176	# if defined (__MWERKS__) & defined (MAC_VERSION)
177	/* Set up interface when using the Metrowerks compiler. This
178	should work for either Macintosh or Windows. */
179	SIOUXSetTitle("\pMrBayes v3.2");
180	SIOUXSettings.fontface = 0; /* plain=0; bold=1 */
181	SIOUXSettings.setupmenus = 0;
182	SIOUXSettings.autocloseonquit = 1;
183	SIOUXSettings.asktosaveonclose = 0;
184	SIOUXSettings.rows = 60;
185	SIOUXSettings.columns = 90;
186	# endif
187
188	# if defined (MPI_ENABLED)
189	# if defined (MAC_VERSION)
190	printf (" Parallel version of\n\n");
191	# else
192	MrBayesPrint (" Parallel version of\n\n");
193	# endif
194	ierror = MPI_Init(&argc, &argv);
195	if (ierror != MPI_SUCCESS)
196	{
197	MrBayesPrint ("%s Problem initializing MPI\n", spacer);
198	exit (1);
199	}
200	ierror = MPI_Comm_size(MPI_COMM_WORLD, &num_procs);
201	if (ierror != MPI_SUCCESS)
202	{
203	MrBayesPrint ("%s Problem getting the number of processors\n", spacer);
204	exit (1);
205	}
206	ierror = MPI_Comm_rank(MPI_COMM_WORLD, &proc_id);
207	if (ierror != MPI_SUCCESS)
208	{
209	MrBayesPrint ("%s Problem getting processors rank\n", spacer);
210	exit (1);
211	}
212	# endif
213
214	#ifdef HAVE_LIBREADLINE
215	rl_attempted_completion_function = readline_completion;
216	#endif
217	/* Set up parameter table. */
218	SetUpParms ();
219
220	/* initialize seed using current time */
221	GetTimeSeed ();
222
223	/* Initialize the variables of the program. */
224	InitializeMrBayes ();
225
226	/* Print the nifty header. */
227	PrintHeader ();
228
229	/* Go to the command line, process any arguments passed to the program
230	and then wait for input. */
231	i = CommandLine (argc, argv);
232
233	# if defined (MPI_ENABLED)
234	MPI_Finalize();
235	# endif
236
237	if (i == ERROR)
238	return (1);
239	else
240	return (0);
241
242	}
243
244
245
246
247
248	int CommandLine (int argc, char **argv)
249
250	{
251
252	int i, message, nProcessedArgs;
253	char cmdStr[CMD_STRING_LENGTH];
254	#ifdef HAVE_LIBREADLINE
255	#ifndef MPI_ENABLED
256	char *cmdStrP;
257	#endif
258	#endif
259	# if defined (MPI_ENABLED)
260	int ierror;
261	# endif
262
263	for(i=0;i<CMD_STRING_LENGTH;i++) cmdStr[i]='\0';
264
265	/* wait for user-input commands */
266	nProcessedArgs = 1; /* first argument is program name and needs not be processed */
267	if (nProcessedArgs < argc)
268	{
269	mode = NONINTERACTIVE; /* when a command is passed into the program, the default is to exit without listening to stdin */
270	autoClose = YES;
271	autoOverwrite = YES;
272	noWarn = YES;
273	quitOnError = YES;
274	}
275	for (;;)
276	{
277	if (nProcessedArgs < argc)
278	{
279	/* we are here only if a command that has been passed
280	into the program remains to be processed */
281	if (nProcessedArgs == 1 && (strcmp(argv[1],"-i") == 0 \|\| strcmp(argv[1],"-I") == 0))
282	{
283	mode = INTERACTIVE;
284	autoClose = NO;
285	autoOverwrite = YES;
286	noWarn = NO;
287	quitOnError = NO;
288	}
289	else
290	sprintf (cmdStr, "Execute %s", argv[nProcessedArgs]);
291	nProcessedArgs++;
292	}
293	else
294	{
295	/* first check if we are in noninteractive mode and quit if so */
296	if (mode == NONINTERACTIVE)
297	{
298	MrBayesPrint ("%s Tasks completed, exiting program because mode is noninteractive\n", spacer);
299	MrBayesPrint ("%s To return control to the command line after completion of file processing, \n", spacer);
300	MrBayesPrint ("%s set mode to interactive with 'mb -i <filename>' (i is for interactive)\n", spacer);
301	MrBayesPrint ("%s or use 'set mode=interactive'\n\n", spacer);
302	return (NO_ERROR);
303	}
304	/* normally, we simply wait at the prompt for a
305	user action */
306	# if defined (MPI_ENABLED)
307	if (proc_id == 0)
308	{
309	/* do not use readline because OpenMPI does not handle it */
310	MrBayesPrint ("MrBayes > ");
311	fflush (stdin);
312	if (fgets (cmdStr, CMD_STRING_LENGTH - 2, stdin) == NULL)
313	{
314	if (feof(stdin))
315	MrBayesPrint ("%s End of File encountered on stdin; quitting\n", spacer);
316	else
317	MrBayesPrint ("%s Could not read command from stdin; quitting\n", spacer);
318	strcpy (cmdStr,"quit;\n");
319	}
320	}
321	ierror = MPI_Bcast (&cmdStr, CMD_STRING_LENGTH, MPI_CHAR, 0, MPI_COMM_WORLD);
322	if (ierror != MPI_SUCCESS)
323	{
324	MrBayesPrint ("%s Problem broadcasting command string\n", spacer);
325	}
326	# else
327	#ifdef HAVE_LIBREADLINE
328	cmdStrP = readline("MrBayes > ");
329	if(cmdStrP!=NULL)
330	{
331	strncpy(cmdStr,cmdStrP,CMD_STRING_LENGTH - 2);
332	if (*cmdStrP)
333	add_history(cmdStrP);
334	free (cmdStrP);
335	}
336	else /* fall through to if (feof(stdin))..*/
337	#else
338	MrBayesPrint ("MrBayes > ");
339	fflush (stdin);
340	if (fgets (cmdStr, CMD_STRING_LENGTH - 2, stdin) == NULL)
341	#endif
342	{
343	if (feof(stdin))
344	MrBayesPrint ("%s End of File encountered on stdin; quitting\n", spacer);
345	else
346	MrBayesPrint ("%s Could not read command from stdin; quitting\n", spacer);
347	strcpy (cmdStr,"quit;\n");
348	}
349	# endif
350	}
351	i = 0;
352	while (cmdStr[i] != '\0' && cmdStr[i] != '\n')
353	i++;
354	cmdStr[i++] = ';';
355	cmdStr[i] = '\0';
356	MrBayesPrint ("\n");
357	if (cmdStr[0] != ';')
358	{
359	/* check that all characters in the string are valid */
360	if (CheckStringValidity (cmdStr) == ERROR)
361	{
362	MrBayesPrint (" Unknown character in command string\n\n");
363	}
364	else
365	{
366	expecting = Expecting(COMMAND);
367	message = ParseCommand (cmdStr);
368
369	if (message == NO_ERROR_QUIT)
370	return (NO_ERROR);
371
372	if (message == ERROR && quitOnError == YES)
373	{
374	MrBayesPrint ("%s Will exit with signal 1 (error) because quitonerror is set to yes\n", spacer);
375	MrBayesPrint ("%s If you want control to be returned to the command line on error,\n", spacer);
376	MrBayesPrint ("%s use 'mb -i <filename>' (i is for interactive) or use 'set quitonerror=no'\n\n", spacer);
377	return (ERROR);
378	}
379
380	# if defined (MPI_ENABLED)
381	ierror = MPI_Barrier (MPI_COMM_WORLD);
382	if (ierror != MPI_SUCCESS)
383	{
384	MrBayesPrint ("%s Problem at command barrier\n", spacer);
385	}
386	# endif
387
388	MrBayesPrint ("\n");
389	}
390	}
391	}
392
393	}
394
395	#ifdef HAVE_LIBREADLINE
396
397	extern char command_generator(const char text, int state);
398
399	char *readline_completion(const char text, int start, int stop) {
400	char matches = (char ) NULL;
401
402	#ifdef COMPLETIONMATCHES
403	if(start == 0)
404	matches = rl_completion_matches (text, command_generator);
405	#endif
406
407	return (matches);
408	}
409	#endif
410
411	void GetTimeSeed (void)
412
413	{
414
415	time_t curTime;
416
417	# if defined (MPI_ENABLED)
418	int ierror;
419
420	if (proc_id == 0)
421	{
422	curTime = time(NULL);
423	globalSeed = (SafeLong)curTime;
424	if (globalSeed < 0)
425	globalSeed = -globalSeed;
426	}
427	ierror = MPI_Bcast(&globalSeed, 1, MPI_LONG, 0, MPI_COMM_WORLD);
428	if (ierror != MPI_SUCCESS)
429	{
430	MrBayesPrint ("%s Problem broadcasting seed\n", spacer);
431	}
432
433	if (proc_id == 0)
434	{
435	/* Note: swapSeed will often be same as globalSeed */
436	curTime = time(NULL);
437	swapSeed = (SafeLong)curTime;
438	if (swapSeed < 0)
439	swapSeed = -swapSeed;
440	}
441	ierror = MPI_Bcast(&swapSeed, 1, MPI_LONG, 0, MPI_COMM_WORLD);
442	if (ierror != MPI_SUCCESS)
443	{
444	MrBayesPrint ("%s Problem broadcasting swap seed\n", spacer);
445	}
446
447	if (proc_id == 0)
448	{
449	/* Note: runIDSeed will often be same as globalSeed */
450	curTime = time(NULL);
451	runIDSeed = (SafeLong)curTime;
452	if (runIDSeed < 0)
453	runIDSeed = -runIDSeed;
454	}
455	ierror = MPI_Bcast(&runIDSeed, 1, MPI_LONG, 0, MPI_COMM_WORLD);
456	if (ierror != MPI_SUCCESS)
457	{
458	MrBayesPrint ("%s Problem broadcasting run ID seed\n", spacer);
459	}
460
461	# else
462	curTime = time(NULL);
463	globalSeed = (SafeLong)curTime;
464	if (globalSeed < 0)
465	globalSeed = -globalSeed;
466
467	/* Note: swapSeed will often be the same as globalSeed */
468	curTime = time(NULL);
469	swapSeed = (SafeLong)curTime;
470	if (swapSeed < 0)
471	swapSeed = -swapSeed;
472
473	/* Note: runIDSeed will often be the same as globalSeed */
474	curTime = time(NULL);
475	runIDSeed = (SafeLong)curTime;
476	if (runIDSeed < 0)
477	runIDSeed = -runIDSeed;
478
479	# endif
480
481	}
482
483
484
485
486
487	int InitializeMrBayes (void)
488
489	{
490	/* this function initializes the program; only call it at the start of execution */
491
492	int i, j, growthFxn[6];
493
494	nBitsInALong = sizeof(SafeLong) * 8; /* global variable: number of bits in a SafeLong */
495	userLevel = STANDARD_USER; /* default user level */
496
497	readWord = NO; /* should we read a word next ? */
498	readComment = NO; /* should we read comments? (used by tree cmd) */
499	fileNameChanged = NO; /* file name changed ? (used by a few commands) */
500	echoMB = YES; /* flag used by Manual to control printing */
501
502	# if defined (MPI_ENABLED)
503	sprintf(manFileName, "commref_mb%sp.txt", VERSION_NUMBER); /* name of command reference file */
504	# else
505	sprintf(manFileName, "commref_mb%s.txt", VERSION_NUMBER); /* name of command reference file */
506	# endif
507
508	for (i=0; i<NUM_ALLOCS; i++) /* set allocated memory to NO */
509	memAllocs[i] = NO;
510	logToFile = NO; /* should screen output be logged to a file */
511	strcpy(logFileName, "log.out"); /* name of the log file */
512	logFileFp = NULL; /* file pointer to log file */
513	replaceLogFile = YES; /* should logfile be replace/appended to */
514	autoClose = NO; /* set default autoclose */
515	autoOverwrite = YES; /* set default autoOverwrite */
516	noWarn = NO; /* set default */
517	quitOnError = NO; /* set default quitOnError */
518	inferAncStates = NO; /* set default inferAncStates */
519	inferSiteOmegas = NO; /* set default inferSiteOmegas */
520	inferSiteRates = NO; /* set default inferSiteRates */
521	inferPosSel = NO; /* set default inferPosSel */
522	inComment = NO; /* not in comment */
523	numComments = 0; /* no comments encountered yet */
524	mode = INTERACTIVE; /* set default mode */
525	numOpenExeFiles = 0; /* no execute files open yet */
526	scientific = YES; /* print to file using scientific format? */
527	precision = 15; /* set default precision */
528	showmovesParams.allavailable = NO; /* do not show all available moves */
529	strcpy(workingDir,""); /* working directory */
530	#if defined (BEAGLE_ENABLED)
531	#if defined (WIN_VERSION)
532	tryToUseBEAGLE = NO; /* try to use the BEAGLE library (NO until SSE code works in Win) */
533	#else
534	tryToUseBEAGLE = NO; /* try to use the BEAGLE library if not Windows (NO untill SSE single prec. works )*/
535	#endif
536	beagleScalingScheme = MB_BEAGLE_SCALE_ALWAYS; /* use BEAGLE dynamic scaling */
537	beagleFlags = BEAGLE_FLAG_PROCESSOR_CPU; /* default to generic CPU */
538	// SSE instructions do not work in Windows environment
539	// beagleFlags \|= BEAGLE_FLAG_VECTOR_SSE; /* default to SSE code */
540	beagleResource = NULL;
541	beagleResourceCount = 0; /* default has no list */
542	beagleInstanceCount = 0; /* no BEAGLE instances */
543	beagleScalingFrequency = 1000;
544	#endif
545	#if defined (THREADS_ENABLED)
546	tryToUseThreads = NO; /* try to use pthread with BEAGLE library */
547	#endif
548
549	/* set the proposal information */
550	SetUpMoveTypes ();
551
552	/* Set up rates for standard amino acid models, in case we need them. */
553	if (SetAARates () == ERROR)
554	return (ERROR);
555
556	/* set up doublet information */
557	doublet[ 0].first = 1; doublet[ 0].second = 1;
558	doublet[ 1].first = 1; doublet[ 1].second = 2;
559	doublet[ 2].first = 1; doublet[ 2].second = 4;
560	doublet[ 3].first = 1; doublet[ 3].second = 8;
561	doublet[ 4].first = 2; doublet[ 4].second = 1;
562	doublet[ 5].first = 2; doublet[ 5].second = 2;
563	doublet[ 6].first = 2; doublet[ 6].second = 4;
564	doublet[ 7].first = 2; doublet[ 7].second = 8;
565	doublet[ 8].first = 4; doublet[ 8].second = 1;
566	doublet[ 9].first = 4; doublet[ 9].second = 2;
567	doublet[10].first = 4; doublet[10].second = 4;
568	doublet[11].first = 4; doublet[11].second = 8;
569	doublet[12].first = 8; doublet[12].second = 1;
570	doublet[13].first = 8; doublet[13].second = 2;
571	doublet[14].first = 8; doublet[14].second = 4;
572	doublet[15].first = 8; doublet[15].second = 8;
573
574	#if defined (FAST_LOG)
575	/* set up log table */
576	for (i=0; i<400; i++)
577	{
578	scalerValue[i] = (CLFlt) ldexp (1.0, 1-i); /* offset 1 needed to deal with scaler == 1.0 */
579	logValue[i] = (CLFlt) log (scalerValue[i]);
580	}
581	#endif
582
583	/* user trees */
584	for (i=0; i<MAX_NUM_USERTREES; i++)
585	userTree[i] = NULL;
586
587	/* parameter values */
588	paramValues = NULL;
589	intValues = NULL;
590
591	/* Prior model settings */
592	defaultModel.dataType = DNA; /* datatype */
593	strcpy(defaultModel.coding, "All"); /* ascertainment bias */
594	strcpy(defaultModel.nucModel, "4by4"); /* nucleotide model */
595	strcpy(defaultModel.nst, "1"); /* number of substitution types */
596	strcpy(defaultModel.aaModelPr, "Fixed"); /* amino acid model prior */
597	for (i=0; i<10; i++)
598	defaultModel.aaModelPrProbs[i] = 0.0;
599	strcpy(defaultModel.aaModel, "Poisson"); /* amino acid model */
600	strcpy(defaultModel.parsModel, "No"); /* do not use parsimony model */
601	strcpy(defaultModel.geneticCode, "Universal"); /* genetic code */
602	strcpy(defaultModel.ploidy, "Diploid"); /* ploidy level */
603	strcpy(defaultModel.omegaVar, "Equal"); /* omega variation */
604	strcpy(defaultModel.ratesModel, "Equal"); /* rates across sites model */
605	defaultModel.numGammaCats = 4; /* number of categories for gamma approximation */
606	strcpy(defaultModel.useGibbs,"No"); /* do not use Gibbs sampling of rate cats by default */
607	defaultModel.gibbsFreq = 100; /* default Gibbs sampling frequency of rate cats*/
608	defaultModel.numBetaCats = 5; /* number of categories for beta approximation */
609	strcpy(defaultModel.covarionModel, "No"); /* use covarion model? (yes/no) */
610	strcpy(defaultModel.augmentData, "No"); /* should data be augmented */
611	strcpy(defaultModel.tRatioPr, "Beta"); /* prior for ti/tv rate ratio */
612	defaultModel.tRatioFix = 1.0;
613	defaultModel.tRatioDir[0] = 1.0;
614	defaultModel.tRatioDir[1] = 1.0;
615	strcpy(defaultModel.revMatPr, "Dirichlet"); /* prior for GTR model (nucleotides) */
616	for (i=0; i<6; i++)
617	{
618	defaultModel.revMatFix[i] = 1.0;
619	defaultModel.revMatDir[i] = 1.0;
620	}
621	defaultModel.revMatSymDir = 1.0; /* default prior for GTR mixed model */
622	strcpy (defaultModel.aaRevMatPr, "Dirichlet"); /* prior for GTR model (proteins) */
623	for (i=0; i<190; i++)
624	{
625	defaultModel.aaRevMatFix[i] = 1.0;
626	defaultModel.aaRevMatDir[i] = 1.0;
627	}
628	strcpy(defaultModel.omegaPr, "Dirichlet"); /* prior for omega */
629	defaultModel.omegaFix = 1.0;
630	defaultModel.omegaDir[0] = 1.0;
631	defaultModel.omegaDir[1] = 1.0;
632	strcpy(defaultModel.ny98omega1pr, "Beta"); /* prior for class 1 omega (Ny98 model) */
633	defaultModel.ny98omega1Fixed = 0.1;
634	defaultModel.ny98omega1Beta[0] = 1.0;
635	defaultModel.ny98omega1Beta[1] = 1.0;
636	strcpy(defaultModel.ny98omega3pr, "Exponential");/* prior for class 3 omega (Ny98 model) */
637	defaultModel.ny98omega3Fixed = 2.0;
638	defaultModel.ny98omega3Uni[0] = 1.0;
639	defaultModel.ny98omega3Uni[1] = 50.0;
640	defaultModel.ny98omega3Exp = 1.0;
641	strcpy(defaultModel.m3omegapr, "Exponential"); /* prior for all three omegas (M3 model) */
642	defaultModel.m3omegaFixed[0] = 0.1;
643	defaultModel.m3omegaFixed[1] = 1.0;
644	defaultModel.m3omegaFixed[2] = 2.0;
645	strcpy(defaultModel.m10betapr, "Uniform"); /* prior for omega variation (M10 model) */
646	strcpy(defaultModel.m10gammapr, "Uniform");
647	defaultModel.m10betaUni[0] = 0.0;
648	defaultModel.m10betaUni[1] = 20.0;
649	defaultModel.m10betaExp = 1.0;
650	defaultModel.m10betaFix[0] = 1.0;
651	defaultModel.m10betaFix[1] = 1.0;
652	defaultModel.m10gammaUni[0] = 0.0;
653	defaultModel.m10gammaUni[1] = 20.0;
654	defaultModel.m10gammaExp = 1.0;
655	defaultModel.m10gammaFix[0] = 1.0;
656	defaultModel.m10gammaFix[1] = 1.0;
657	defaultModel.numM10GammaCats = 4;
658	defaultModel.numM10BetaCats = 4;
659	strcpy(defaultModel.codonCatFreqPr, "Dirichlet");/* prior for selection cat frequencies */
660	defaultModel.codonCatFreqFix[0] = 1.0/3.0;
661	defaultModel.codonCatFreqFix[1] = 1.0/3.0;
662	defaultModel.codonCatFreqFix[2] = 1.0/3.0;
663	defaultModel.codonCatDir[0] = 1.0;
664	defaultModel.codonCatDir[1] = 1.0;
665	defaultModel.codonCatDir[2] = 1.0;
666	strcpy(defaultModel.stateFreqPr, "Dirichlet"); /* prior for character state frequencies */
667	strcpy(defaultModel.stateFreqsFixType, "Equal");
668	for (i=0; i<200; i++)
669	{
670	defaultModel.stateFreqsFix[i] = 0.0;
671	defaultModel.stateFreqsDir[i] = 1.0;
672	}
673	defaultModel.numDirParams = 0;
674	strcpy(defaultModel.shapePr, "Uniform"); /* prior for gamma shape parameter */
675	defaultModel.shapeFix = 0.5;
676	defaultModel.shapeUni[0] = MIN_SHAPE_PARAM;
677	defaultModel.shapeUni[1] = MAX_SHAPE_PARAM;
678	defaultModel.shapeExp = 2.0;
679	strcpy(defaultModel.pInvarPr, "Uniform"); /* prior for proportion of invariable sites */
680	defaultModel.pInvarFix = 0.1;
681	defaultModel.pInvarUni[0] = 0.0;
682	defaultModel.pInvarUni[1] = 1.0;
683	strcpy(defaultModel.adGammaCorPr, "Uniform"); /* prior for correlation param of adGamma model */
684	defaultModel.corrFix = 0.0;
685	defaultModel.corrUni[0] = -1.0;
686	defaultModel.corrUni[1] = 1.0;
687	strcpy(defaultModel.covSwitchPr, "Uniform"); /* prior for switching rates of covarion model */
688	defaultModel.covswitchFix[0] = 1.0;
689	defaultModel.covswitchFix[1] = 1.0;
690	defaultModel.covswitchUni[0] = 0.0;
691	defaultModel.covswitchUni[1] = 100.0;
692	defaultModel.covswitchExp = 1.0;
693	strcpy(defaultModel.symPiPr, "Fixed"); /* prior for pi when unidentifiable states used */
694	defaultModel.symBetaFix = -1.0;
695	defaultModel.symBetaUni[0] = 0.0;
696	defaultModel.symBetaUni[1] = 20.0;
697	defaultModel.symBetaExp = 2;
698	strcpy(defaultModel.brownCorPr, "Fixed"); /* prior on correlation of brownian model */
699	defaultModel.brownCorrFix = 0.0;
700	defaultModel.brownCorrUni[0] = -1.0;
701	defaultModel.brownCorrUni[1] = 1.0;
702	strcpy(defaultModel.brownScalesPr, "Gammamean"); /* prior on scales of brownian model */
703	defaultModel.brownScalesFix = 10.0;
704	defaultModel.brownScalesUni[0] = 0.0;
705	defaultModel.brownScalesUni[1] = 100.0;
706	defaultModel.brownScalesGamma[0] = 1.0;
707	defaultModel.brownScalesGamma[1] = 10.0;
708	defaultModel.brownScalesGammaMean = 10.0;
709	strcpy(defaultModel.topologyPr, "Uniform"); /* prior for tree topology */
710	defaultModel.topologyFix = -1; /* user tree index to use for fixed topology */
711	defaultModel.activeConstraints = NULL; /* which constraints are active */
712	strcpy(defaultModel.brlensPr, "Unconstrained"); /* prior on branch lengths */
713	defaultModel.brlensFix = -1; /* user tree index to use for fixed brlens */
714	defaultModel.brlensUni[0] = BRLENS_MIN;
715	defaultModel.brlensUni[1] = 10.0;
716	defaultModel.brlensExp = 10.0;
717	strcpy(defaultModel.unconstrainedPr, "Exponential");/* prior on branches if unconstrained */
718	strcpy(defaultModel.clockPr, "Uniform"); /* prior on branch lengths if clock enforced */
719	strcpy(defaultModel.treeAgePr, "Exponential"); /* prior on tree age */
720	defaultModel.treeAgeGamma[0] = 1.0;
721	defaultModel.treeAgeGamma[1] = 1.0;
722	defaultModel.treeAgeExp = 1.0;
723	defaultModel.treeAgeFix = 1.0;
724	strcpy(defaultModel.clockRatePr, "Fixed"); /* prior on base subst. rate for clock trees */
725	defaultModel.clockRateNormal[0] = 1.0;
726	defaultModel.clockRateNormal[1] = 1.0;
727	defaultModel.clockRateLognormal[0] = 0.0; /* mean 0.0 on log scale corresponds to mean rate 1.0 */
728	defaultModel.clockRateLognormal[1] = 0.7; /* double or half the rate in one standard deviation */
729	defaultModel.clockRateGamma[0] = 1.0;
730	defaultModel.clockRateGamma[1] = 1.0;
731	defaultModel.clockRateExp = 1.0;
732	defaultModel.clockRateFix = 1.0;
733	strcpy(defaultModel.speciationPr, "Exponential"); /* prior on speciation rate */
734	defaultModel.speciationFix = 1.0;
735	defaultModel.speciationUni[0] = 0.0;
736	defaultModel.speciationUni[1] = 10.0;
737	defaultModel.speciationExp = 1.0;
738	strcpy(defaultModel.extinctionPr, "Beta"); /* prior on extinction rate */
739	defaultModel.extinctionFix = 1.0;
740	defaultModel.extinctionBeta[0] = 1;
741	defaultModel.extinctionBeta[1] = 1;
742	strcpy(defaultModel.sampleStrat, "Random"); /* taxon sampling strategy */
743	defaultModel.sampleProb = 1.0; /* taxon sampling fraction */
744	strcpy(defaultModel.popSizePr, "Lognormal"); /* prior on coalescence population size */
745	defaultModel.popSizeFix = 10.0;
746	defaultModel.popSizeUni[0] = 1.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] = 2.3; /* factor 10 in one standard deviation */
752	defaultModel.popSizeGamma[0] = 100.0;
753	defaultModel.popSizeGamma[1] = 1000.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 = 10.0;
771	defaultModel.tk02varFix = 0.1;
772	defaultModel.tk02varUni[0] = 0.0;
773	defaultModel.tk02varUni[1] = 0.5;
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.ratePr, "Fixed"); /* prior on rate for a partition */
780	defaultModel.ratePrDir = 1.0;
781
782	defaultModel.nStates = 4; /* number of states for partition */
783
784	/* Report format settings */
785	strcpy(defaultModel.tratioFormat, "Ratio"); /* default format for tratio */
786	strcpy(defaultModel.revmatFormat, "Dirichlet"); /* default format for revmat */
787	strcpy(defaultModel.ratemultFormat, "Scaled"); /* default format for ratemult */
788	strcpy(defaultModel.treeFormat, "Brlens"); /* default format for trees */
789	strcpy(defaultModel.inferAncStates, "No"); /* do not infer ancestral states */
790	strcpy(defaultModel.inferPosSel, "No"); /* do not infer positive selection */
791	strcpy(defaultModel.inferSiteOmegas, "No"); /* do not infer site omega vals */
792	strcpy(defaultModel.inferSiteRates, "No"); /* do not infer site rates */
793
794	/* Allocate and initialize model indicator parameter names */
795	modelIndicatorParams = (char *) SafeCalloc (3, sizeof (char ));
796	modelIndicatorParams[0] = "aamodel";
797	modelIndicatorParams[1] = "gtrsubmodel";
798	modelIndicatorParams[2] = "";
799
800	/* Aamodel */
801	modelElementNames = (char *) SafeCalloc (3, sizeof (char ));
802	modelElementNames[0] = (char *) SafeCalloc (11, sizeof (char ));
803	modelElementNames[0][0] = "Poisson";
804	modelElementNames[0][1] = "Jones";
805	modelElementNames[0][2] = "Dayhoff";
806	modelElementNames[0][3] = "Mtrev";
807	modelElementNames[0][4] = "Mtmam";
808	modelElementNames[0][5] = "Wag";
809	modelElementNames[0][6] = "Rtrev";
810	modelElementNames[0][7] = "Cprev";
811	modelElementNames[0][8] = "Vt";
812	modelElementNames[0][9] = "Blosum";
813	modelElementNames[0][10] = "";
814
815	/* Gtrsubmodel */
816	modelElementNames[1] = (char *) SafeCalloc (204, sizeof (char ));
817	for (i=0; i<203; i++)
818	{
819	modelElementNames[1][i] = (char *) SafeCalloc (7, sizeof (char));
820	FromIndexToGrowthFxn(i, growthFxn);
821	for (j=0; j<6; j++)
822	modelElementNames[1][i][j] = '1' + growthFxn[j];
823	modelElementNames[1][i][j] = '\0';
824	}
825	modelElementNames[1][203] = "";
826
827
828	/* Termination */
829	modelElementNames[2] = (char *) SafeCalloc (1, sizeof(char ));
830	modelElementNames[2][0] = "";
831
832	/* initialize user trees */
833	for (i=0; i<MAX_NUM_USERTREES; i++)
834	userTree[i] = NULL;
835	numUserTrees = 0;
836
837	/* Reset translate table */
838	ResetTranslateTable();
839
840	/* finally reset everything dependent on a matrix being defined */
841	return (ReinitializeMrBayes ());
842
843	}
844
845
846
847	unsigned FindMaxRevision ( unsigned amount, ...)
848	{
849	const char* cur;
850	char tmp[20];
851	unsigned val,i,max;
852
853	va_list vl;
854	va_start(vl,amount);
855	max=0;
856	for (i=0;i<amount;i++)
857	{
858	cur=va_arg(vl,const char*);
859	sscanf(cur,"%s %d",tmp,&val);
860	max=(max>val)?max:val;
861	}
862	va_end(vl);
863	return max;
864	}
865
866
867
868	void PrintHeader (void)
869
870	{
871	char arch[4];
872	#ifndef RELEASE
873	unsigned rev=FindMaxRevision ( 13, svnRevisionBayesC,svnRevisionBestC,svnRevisionCommandC,svnRevisionMbC,svnRevisionMbbeagleC,svnRevisionMbmathC,svnRevisionMcmcC,svnRevisionModelC,svnRevisionPlotC,svnRevisionSumpC,svnRevisionSumtC,svnRevisionTreeC,svnRevisionUtilsC);
874	#endif
875
876	strcpy(arch,(sizeof(void*)==4)?"x86":"x64");
877
878	# if defined (MAC_VERSION)
879
880	# if !defined (MPI_ENABLED)
881	printf ("\n\n");
882	# endif
883	#ifdef RELEASE
884	printf (" MrBayes v%s %s\n\n", VERSION_NUMBER,arch);
885	#else
886	printf (" MrBayes v%s(r%d) %s\n\n", VERSION_NUMBER,rev,arch);
887	#endif
888	printf (" (Bayesian Analysis of Phylogeny)\n\n");
889	# if defined (MPI_ENABLED)
890	printf (" (Parallel version)\n");
891	printf (" (%d processors available)\n\n", num_procs);
892	# endif
893	printf (" Distributed under the GNU General Public License\n\n\n");
894
895	# if defined (MPI_ENABLED)
896	if (proc_id == 0)
897	{
898	printf (" Node number %d\n\n", proc_id + 1);
899	printf (" Type \"help\" or \"help <command>\" for information\n");
900	printf (" on the commands that are available.\n\n\n");
901	}
902	else
903	{
904	printf (" Remote node number %d\n\n", proc_id + 1);
905	printf (" All information is printed to node 1\n\n\n");
906	}
907	# else
908	printf (" Type \"help\" or \"help <command>\" for information\n");
909	printf (" on the commands that are available.\n\n");
910	printf (" Type \"about\" for authorship and general\n");
911	printf (" information about the program\n\n\n");
912	# endif
913
914	# else
915
916	# if !defined (MPI_ENABLED)
917	MrBayesPrint ("\n\n");
918	# endif
919	#ifdef RELEASE
920	MrBayesPrint (" MrBayes v%s %s\n\n", VERSION_NUMBER,arch);
921	#else
922	MrBayesPrint (" MrBayes v%s(r%d) %s\n\n", VERSION_NUMBER,rev,arch);
923	#endif
924	MrBayesPrint (" (Bayesian Analysis of Phylogeny)\n\n");
925	# if defined (MPI_ENABLED)
926	MrBayesPrint (" (Parallel version)\n");
927	MrBayesPrint (" (%d processors available)\n\n", num_procs);
928	# endif
929	MrBayesPrint (" Distributed under the GNU General Public License\n\n\n");
930	MrBayesPrint (" Type \"help\" or \"help <command>\" for information\n");
931	MrBayesPrint (" on the commands that are available.\n\n");
932	MrBayesPrint (" Type \"about\" for authorship and general\n");
933	MrBayesPrint (" information about the program.\n\n\n");
934	# endif
935
936	}
937
938
939
940
941
942	int ReinitializeMrBayes (void)
943
944	{
945
946	/* this function resets everything dependent on a matrix */
947
948	int i;
949
950	/* reinitialize indentation */
951	strcpy (spacer, ""); /* holds blanks for indentation */
952
953	/* reset all taxa flags */
954	ResetTaxaFlags();
955
956	/* reset all characters flags */
957	ResetCharacterFlags();
958
959	/* chain parameters */
960	chainParams.numGen = 1000000; /* number of MCMC cycles */
961	chainParams.sampleFreq = 500; /* frequency to sample chain */
962	chainParams.printFreq = 500; /* frequency to print chain */
963	chainParams.swapFreq = 1; /* frequency of attempting swap of states */
964	chainParams.numSwaps = 1; /* number of swaps to try each time */
965	chainParams.isSS = NO;
966	chainParams.startFromPriorSS = NO;
967	chainParams.numStepsSS = 50;
968	chainParams.burninSS = -1;
969	chainParams.alphaSS = 0.4;
970	chainParams.backupCheckSS = 0;
971	chainParams.mcmcDiagn = YES; /* write MCMC diagnostics to file ? */
972	chainParams.diagnFreq = 5000; /* diagnostics frequency */
973	chainParams.minPartFreq = 0.10; /* min partition frequency for diagnostics */
974	chainParams.allChains = NO; /* calculate diagnostics for all chains ? */
975	chainParams.allComps = NO; /* do not calc diagn for all run comparisons */
976	chainParams.relativeBurnin = YES; /* use relative burnin? */
977	chainParams.burninFraction = 0.25; /* default burnin fraction */
978	chainParams.stopRule = NO; /* should stopping rule be used? */
979	chainParams.stopVal = 0.05; /* convergence diagnostic value to reach */
980	chainParams.numRuns = 2; /* number of runs */
981	chainParams.numChains = 4; /* number of chains */
982	chainParams.chainTemp = 0.1; /* chain temperature */
983	chainParams.redirect = NO; /* should printf be to stdout */
984	strcpy(chainParams.chainFileName, "temp.out"); /* chain file name for output */
985	chainParams.chainBurnIn = 0; /* chain burn in length */
986	chainParams.numStartPerts = 0; /* number of perturbations to starting tree */
987	strcpy(chainParams.startTree, "Current"); /* starting tree for chain (random/current) */
988	strcpy(chainParams.startParams, "Current"); /* starting params for chain (reset/current) */
989	chainParams.saveBrlens = YES; /* should branch lengths be saved */
990	chainParams.weightScheme[0] = 0.0; /* percent chars to decrease in weight */
991	chainParams.weightScheme[1] = 0.0; /* percent chars to increase in weight */
992	chainParams.weightScheme[2] = 1.0; /* weight increment */
993	chainParams.calcPbf = NO; /* should we calculate the pseudo-BF? */
994	chainParams.pbfInitBurnin = 100000; /* initial burnin for pseudo BF */
995	chainParams.pbfSampleFreq = 10; /* sample frequency for pseudo BF */
996	chainParams.pbfSampleTime = 2000; /* how many cycles to calcualate site prob. */
997	chainParams.pbfSampleBurnin = 2000; /* burnin period for each site for pseudo BF */
998	chainParams.userDefinedTemps = NO; /* should we use the users temperatures? */
999	for (i=0; i<MAX_CHAINS; i++)
1000	chainParams.userTemps[i] = 1.0; /* user-defined chain temperatures */
1001	chainParams.swapAdjacentOnly = NO; /* swap only adjacent temperatures */
1002	chainParams.printMax = 8; /* maximum number of chains to print to screen */
1003	chainParams.printAll = YES; /* whether to print heated chains */
1004	chainParams.treeList = NULL; /* vector of tree lists for saving trees */
1005	chainParams.saveTrees = NO; /* save tree samples for later removal? */
1006	chainParams.tFilePos = NULL; /* position for reading trees for removal */
1007	chainParams.runWithData = YES; /* whether to run with data */
1008	chainParams.orderTaxa = NO; /* should taxa be ordered in output trees? */
1009	chainParams.append = NO; /* append to previous analysis? */
1010	chainParams.autotune = YES; /* autotune? */
1011	chainParams.tuneFreq = 100; /* autotuning frequency */
1012	chainParams.checkPoint = YES; /* should we checkpoint the run? */
1013	chainParams.checkFreq = 100000; /* check-pointing frequency */
1014	chainParams.diagnStat = AVGSTDDEV; /* mcmc diagnostic to use */
1015
1016	/* sumt parameters */
1017	strcpy(sumtParams.sumtFileName, "temp.t"); /* input name for sumt command */
1018	strcpy(sumtParams.sumtConType, "Halfcompat"); /* type of consensus tree output */
1019	sumtParams.calcTreeprobs = YES; /* should individual tree probs be calculated */
1020	sumtParams.showSumtTrees = NO; /* should the individual tree probs be shown */
1021	sumtParams.printBrlensToFile = NO; /* should brlens be printed to file */
1022	sumtParams.brlensFreqDisplay = 0.50; /* threshold for printing brlens to file */
1023	sumtParams.numTrees = 1; /* number of trees to summarize */
1024	sumtParams.numRuns = 2; /* number of analyses to summarize */
1025	sumtParams.orderTaxa = YES; /* order taxa in trees ? */
1026	sumtParams.minPartFreq = 0.10; /* minimum part. freq. for overall diagnostics */
1027	sumtParams.table = YES; /* display table of part. freq.? */
1028	sumtParams.summary = YES; /* display overall diagnostics? */
1029	sumtParams.showConsensus = YES; /* display consensus tree(s)? */
1030	sumtParams.consensusFormat = FIGTREE; /* format of consensus tree */
1031	strcpy (sumtParams.sumtOutfile, "temp.t.stat"); /* output name for sumt command */
1032	sumtParams.HPD = YES; /* use Highest Posterior Density? */
1033
1034	/* sump parameters */
1035	strcpy(sumpParams.sumpFileName, "temp.p"); /* input name for sump command */
1036	sumpParams.numRuns = 2; /* number of analyses to summarize */
1037	sumpParams.HPD = YES; /* use Highest Posterior Density? */
1038	sumpParams.minProb = 0.05; /* min. prob. of models to include in summary */
1039
1040	/* sumss parameters */
1041	sumssParams.numRuns= 2; /* number of independent analyses to summarize */
1042	sumssParams.allRuns = YES; /* should data for all runs be printed (yes/no)? */
1043	sumssParams.stepToPlot = 0; /* Which step to plot in the step plot, 0 means burnin */
1044	sumssParams.askForMorePlots = YES; /* Should user be asked to plot for different discardfraction (y/n)? */
1045	sumssParams.smoothing = 0; /* An integer indicating number of neighbors to average over when dooing smoothing of curvs on plots */
1046	sumssParams.discardFraction = 0.8; /* Proportion of samples discarded when ploting step plot.*/
1047
1048	/* comparetree parameters */
1049	strcpy(comptreeParams.comptFileName1, "temp.t"); /* input name for comparetree command */
1050	strcpy(comptreeParams.comptFileName2, "temp.t"); /* input name for comparetree command */
1051	strcpy(comptreeParams.comptOutfile, "temp.comp");/* input name for comparetree command */
1052	comptreeParams.minPartFreq = 0.0; /* minimum frequency of partitions to include */
1053
1054	/* plot parameters */
1055	strcpy(plotParams.plotFileName, "temp.p"); /* input name for plot command */
1056	strcpy(plotParams.parameter, "lnL"); /* plotted parameter plot command */
1057	strcpy(plotParams.match, "Perfect"); /* matching for plot command */
1058
1059	return (NO_ERROR);
1060
1061	}
1062
1063
1064
1065
1066
1067	int DoQuit (void)
1068
1069	{
1070
1071	int i;
1072	char tempName[100];
1073
1074	/* free information for model and matrix */
1075	FreeModel ();
1076	FreeMatrix ();
1077
1078	SafeFclose (&logFileFp);
1079	logToFile = NO;
1080
1081	/* check to see if any memory has not been freed */
1082	for (i=0; i<NUM_ALLOCS; i++)
1083	{
1084	if (memAllocs[i] == YES)
1085	{
1086	MrBayesPrint (" WARNING: Memory (%d) has not been freed\n", i);
1087	if (mode == INTERACTIVE && quitOnError == NO)
1088	{
1089	MrBayesPrint ("%s Hit return key to continue ", spacer);
1090	fflush (stdin);
1091	if( fgets (tempName, 100, stdin) == NULL )
1092	{
1093	printf("Error in function: %s at line: %d in file: %s", __FUNCTION__, __LINE__, __FILE__);
1094	}
1095	}
1096	}
1097	}
1098
1099	/* free modelIndicatorParams and modelElementNames */
1100	for (i=0; i<203; i++)
1101	free (modelElementNames[1][i]);
1102	for (i=0; i<3; i++)
1103	free (modelElementNames[i]);
1104	free (modelElementNames);
1105	free (modelIndicatorParams);
1106
1107	MrBayesPrint (" Quitting program\n\n");
1108
1109	/* If we quit while reading a mrbayes block, then we need to make certain
1110	that we return a NO_ERROR_QUIT so we can break out of DoExecute cleanly,
1111	and dealloc "s" there. */
1112	if (inMrbayesBlock == YES)
1113	{
1114	inMrbayesBlock = NO;
1115	return (NO_ERROR_QUIT);
1116	}
1117
1118	return (NO_ERROR);
1119
1120	}
1121
1122
1123
1124
1125
1126	void SetCode (int part)
1127
1128	{
1129
1130	int i, s, s1, s2, s3, ns;
1131
1132	modelParams[part].codon[ 0] = 12; /* AAA Lys */
1133	modelParams[part].codon[ 1] = 3; /* AAC Asn */
1134	modelParams[part].codon[ 2] = 12; /* AAG Lys */
1135	modelParams[part].codon[ 3] = 3; /* AAT Asn */
1136	modelParams[part].codon[ 4] = 17; /* ACA Thr */
1137	modelParams[part].codon[ 5] = 17; /* ACC Thr */
1138	modelParams[part].codon[ 6] = 17; /* ACG Thr */
1139	modelParams[part].codon[ 7] = 17; /* ACT Thr */
1140	modelParams[part].codon[ 8] = 2; /* AGA Arg */
1141	modelParams[part].codon[ 9] = 16; /* AGC Ser */
1142	modelParams[part].codon[10] = 2; /* AGG Arg */
1143	modelParams[part].codon[11] = 16; /* AGT Ser */
1144	modelParams[part].codon[12] = 10; /* ATA Ile */
1145	modelParams[part].codon[13] = 10; /* ATC Ile */
1146	modelParams[part].codon[14] = 13; /* ATG Met */
1147	modelParams[part].codon[15] = 10; /* ATT Ile */
1148	modelParams[part].codon[16] = 6; /* CAA Gln */
1149	modelParams[part].codon[17] = 9; /* CAC His */
1150	modelParams[part].codon[18] = 6; /* CAG Gln */
1151	modelParams[part].codon[19] = 9; /* CAT His */
1152	modelParams[part].codon[20] = 15; /* CCA Pro */
1153	modelParams[part].codon[21] = 15; /* CCC Pro */
1154	modelParams[part].codon[22] = 15; /* CCG Pro */
1155	modelParams[part].codon[23] = 15; /* CCT Pro */
1156	modelParams[part].codon[24] = 2; /* CGA Arg */
1157	modelParams[part].codon[25] = 2; /* CGC Arg */
1158	modelParams[part].codon[26] = 2; /* CGG Arg */
1159	modelParams[part].codon[27] = 2; /* CGT Arg */
1160	modelParams[part].codon[28] = 11; /* CTA Leu */
1161	modelParams[part].codon[29] = 11; /* CTC Leu */
1162	modelParams[part].codon[30] = 11; /* CTG Leu */
1163	modelParams[part].codon[31] = 11; /* CTT Leu */
1164	modelParams[part].codon[32] = 7; /* GAA Glu */
1165	modelParams[part].codon[33] = 4; /* GAC Asp */
1166	modelParams[part].codon[34] = 7; /* GAG Glu */
1167	modelParams[part].codon[35] = 4; /* GAT Asp */
1168	modelParams[part].codon[36] = 1; /* GCA Ala */
1169	modelParams[part].codon[37] = 1; /* GCC Ala */
1170	modelParams[part].codon[38] = 1; /* GCG Ala */
1171	modelParams[part].codon[39] = 1; /* GCT Ala */
1172	modelParams[part].codon[40] = 8; /* GGA Gly */
1173	modelParams[part].codon[41] = 8; /* GGC Gly */
1174	modelParams[part].codon[42] = 8; /* GGG Gly */
1175	modelParams[part].codon[43] = 8; /* GGT Gly */
1176	modelParams[part].codon[44] = 20; /* GTA Val */
1177	modelParams[part].codon[45] = 20; /* GTC Val */
1178	modelParams[part].codon[46] = 20; /* GTG Val */
1179	modelParams[part].codon[47] = 20; /* GTT Val */
1180	modelParams[part].codon[48] = 21; /* TAA Stop*/
1181	modelParams[part].codon[49] = 19; /* TAC Tyr */
1182	modelParams[part].codon[50] = 21; /* TAG Stop*/
1183	modelParams[part].codon[51] = 19; /* TAT Tyr */
1184	modelParams[part].codon[52] = 16; /* TCA Ser */
1185	modelParams[part].codon[53] = 16; /* TCC Ser */
1186	modelParams[part].codon[54] = 16; /* TCG Ser */
1187	modelParams[part].codon[55] = 16; /* TCT Ser */
1188	modelParams[part].codon[56] = 21; /* TGA Stop*/
1189	modelParams[part].codon[57] = 5; /* TGC Cys */
1190	modelParams[part].codon[58] = 18; /* TGG Trp */
1191	modelParams[part].codon[59] = 5; /* TGT Cys */
1192	modelParams[part].codon[60] = 11; /* TTA Leu */
1193	modelParams[part].codon[61] = 14; /* TTC Phe */
1194	modelParams[part].codon[62] = 11; /* TTG Leu */
1195	modelParams[part].codon[63] = 14; /* TTT Phe */
1196
1197	if (!strcmp(modelParams[part].geneticCode, "Vertmt"))
1198	{
1199	/* UGA: Ter -> Trp
1200	AUA: Ile -> Met
1201	AGA: Arg -> Ter
1202	AGG: Arg -> Ter */
1203	modelParams[part].codon[ 8] = 21; /* AGA Stop */
1204	modelParams[part].codon[10] = 21; /* AGG Stop */
1205	modelParams[part].codon[12] = 13; /* ATA Met */
1206	modelParams[part].codon[56] = 18; /* TGA Trp */
1207	}
1208	else if (!strcmp(modelParams[part].geneticCode, "Mycoplasma"))
1209	{
1210	/* UGA: Ter -> Trp */
1211	modelParams[part].codon[56] = 18; /* TGA Trp */
1212	}
1213	else if (!strcmp(modelParams[part].geneticCode, "Yeast"))
1214	{
1215	/* UGA: Ter -> Trp
1216	AUA: Ile -> Met
1217	CUA: Leu -> Thr
1218	CUC: Leu -> Thr
1219	CUG: Leu -> Thr
1220	CUU: Leu -> Thr */
1221	modelParams[part].codon[12] = 13; /* ATA Met */
1222	modelParams[part].codon[28] = 17; /* CTA Thr */
1223	modelParams[part].codon[29] = 17; /* CTC Thr */
1224	modelParams[part].codon[30] = 17; /* CTG Thr */
1225	modelParams[part].codon[31] = 17; /* CTT Thr */
1226	modelParams[part].codon[56] = 18; /* TGA Trp */
1227	}
1228	else if (!strcmp(modelParams[part].geneticCode, "Ciliates"))
1229	{
1230	/* UAA: Ter -> Gln
1231	UAG: Ter -> Gln */
1232	modelParams[part].codon[48] = 6; /* TAA Gln */
1233	modelParams[part].codon[50] = 6; /* TAG Gln */
1234	}
1235	else if (!strcmp(modelParams[part].geneticCode, "Metmt"))
1236	{
1237	/* UGA: Ter -> Trp
1238	AUA: Ile -> Met
1239	AGA: Arg -> Ser
1240	AGG: Arg -> Ser */
1241	modelParams[part].codon[ 8] = 16; /* AGA Ser */
1242	modelParams[part].codon[10] = 16; /* AGG Ser */
1243	modelParams[part].codon[12] = 13; /* ATA Met */
1244	modelParams[part].codon[56] = 18; /* TGA Trp */
1245	}
1246	else
1247	{
1248
1249	}
1250
1251	ns = 0;
1252	for (i=0; i<64; i++)
1253	{
1254	if (modelParams[part].codon[i] != 21)
1255	ns++;
1256	}
1257	/* printf ("ns = %d\n", ns); */
1258
1259	s = 0;
1260	for (s1=0; s1<4; s1++)
1261	{
1262	for (s2=0; s2<4; s2++)
1263	{
1264	for (s3=0; s3<4; s3++)
1265	{
1266	if (modelParams[part].codon[s116 + s24 + s3] != 21)
1267	{
1268	modelParams[part].codonNucs[s][0] = s1;
1269	modelParams[part].codonNucs[s][1] = s2;
1270	modelParams[part].codonNucs[s][2] = s3;
1271	modelParams[part].codonAAs[s] = modelParams[part].codon[s116 + s24 + s3];
1272	s++;
1273	}
1274	}
1275	}
1276	}
1277
1278	}
1279

Note: See TracBrowser for help on using the repository browser.

Context Navigation

source: branches/items/GDE/MrBAYES/mrbayes_3.2.1/bayes.c

Download in other formats: