Context Navigation

plot.c

Visit:

Last change on this file was 10416, checked in by aboeckma, 11 years ago
Added mr bayes (no menu yet)
File size: 9.8 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 <ctype.h>
43	#include "mb.h"
44	#include "globals.h"
45	#include "command.h"
46	#include "bayes.h"
47	#include "plot.h"
48	#include "sump.h"
49	#include "utils.h"
50	#if defined(__MWERKS__)
51	#include "SIOUX.h"
52	#endif
53
54	const char* const svnRevisionPlotC="$Rev: 343 $"; /* Revision keyword which is expended/updated by svn on each commit/update*/
55
56	/* local (to this file) */
57	int foundCurly;
58	char plotTokenP;/ plotToken[CMD_STRING_LENGTH];*/
59
60
61
62
63
64
65	int DoPlot (void)
66
67	{
68
69
70	int i, n, nHeaders, burnin, len, longestHeader, whichIsX, whichIsY, numPlotted;
71	char temp[100], **headerNames = NULL;
72	SumpFileInfo fileInfo;
73	ParameterSample *parameterSamples;
74
75	# if defined (MPI_ENABLED)
76	if (proc_id != 0)
77	return NO_ERROR;
78	# endif
79
80	/* initialize values */
81	headerNames = NULL;
82	nHeaders = 0;
83	parameterSamples = NULL;
84
85	/* tell user we are ready to go */
86	MrBayesPrint ("%s Plotting parameters in file %s ...\n", spacer, plotParams.plotFileName);
87
88	/* examine plot file */
89	if (ExamineSumpFile (plotParams.plotFileName, &fileInfo, &headerNames, &nHeaders) == ERROR)
90	return ERROR;
91
92	/* Calculate burn in */
93	burnin = fileInfo.firstParamLine - fileInfo.headerLine - 1;
94
95	/* tell the user that everything is fine */
96	MrBayesPrint ("%s Found %d parameter lines in file \"%s\"\n", spacer, fileInfo.numRows + burnin, plotParams.plotFileName);
97	if (burnin > 0)
98	MrBayesPrint ("%s Of the %d lines, %d of them will be summarized (starting at line %d)\n", spacer, fileInfo.numRows+burnin, fileInfo.numRows, fileInfo.firstParamLine);
99	else
100	MrBayesPrint ("%s All %d lines will be summarized (starting at line %d)\n", spacer, fileInfo.numRows, fileInfo.firstParamLine);
101	MrBayesPrint ("%s (Only the last set of lines will be read, in case multiple\n", spacer);
102	MrBayesPrint ("%s parameter blocks are present in the same file.)\n", spacer);
103
104	/* allocate space to hold parameter information */
105	if (AllocateParameterSamples (&parameterSamples, 1, fileInfo.numRows, fileInfo.numColumns) == ERROR)
106	goto errorExit;
107
108	/* Now we read the file for real. First, rewind file pointer to beginning of file... */
109	if (ReadParamSamples (plotParams.plotFileName, &fileInfo, parameterSamples, 0) == ERROR)
110	goto errorExit;
111
112	/* get length of longest header */
113	longestHeader = 9; /* 9 is the length of the word "parameter" (for printing table) */
114	for (i=0; i<nHeaders; i++)
115	{
116	len = (int) strlen(headerNames[i]);
117	if (len > longestHeader)
118	longestHeader = len;
119	}
120
121	/* print x-y plot of parameter vs. generation */
122	whichIsX = -1;
123	for (i=0; i<nHeaders; i++)
124	{
125	if (IsSame (headerNames[i], "Gen") == SAME)
126	whichIsX = i;
127	}
128
129	if (whichIsX < 0)
130	{
131	MrBayesPrint ("%s Could not find a column labelled \"Gen\" \n", spacer);
132	goto errorExit;
133	}
134
135	numPlotted = 0;
136	for (n=0; n<nHeaders; n++)
137	{
138	strcpy (temp, headerNames[n]);
139	whichIsY = -1;
140	if (!strcmp(plotParams.match, "Perfect"))
141	{
142	if (IsSame (temp, plotParams.parameter) == SAME)
143	whichIsY = n;
144	}
145	else if (!strcmp(plotParams.match, "All"))
146	{
147	whichIsY = n;
148	}
149	else
150	{
151	if (IsSame (temp, plotParams.parameter) == CONSISTENT_WITH)
152	whichIsY = n;
153	}
154
155	if (whichIsY >= 0 && whichIsX != whichIsY)
156	{
157	MrBayesPrint ("\n%s Rough trace plot of parameter %s:\n", spacer, headerNames[whichIsY]);
158	if (PrintPlot (parameterSamples[whichIsX].values[0], parameterSamples[whichIsY].values[0], fileInfo.numRows) == ERROR)
159	goto errorExit;
160	numPlotted++;
161	}
162	}
163
164	if (numPlotted == 0)
165	{
166	MrBayesPrint ("%s Did not find any parameters matching \"%s\" to plot\n", spacer, plotParams.parameter);
167	}
168
169	/* free memory */
170	for (i=0; i<nHeaders; i++)
171	free (headerNames[i]);
172	free(headerNames);
173	FreeParameterSamples(parameterSamples);
174
175	expecting = Expecting(COMMAND);
176
177	return (NO_ERROR);
178
179	errorExit:
180
181	/* free memory */
182	for (i=0; i<nHeaders; i++)
183	free (headerNames[i]);
184	free(headerNames);
185	FreeParameterSamples(parameterSamples);
186
187	expecting = Expecting(COMMAND);
188
189	return (ERROR);
190	}
191
192
193
194
195
196	int DoPlotParm (char parmName, char tkn)
197
198	{
199
200	int tempI;
201	MrBFlt tempD;
202	char tempStr[100];
203
204	if (defMatrix == NO)
205	{
206	MrBayesPrint ("%s A matrix must be specified before sumt can be used\n", spacer);
207	return (ERROR);
208	}
209
210	if (expecting == Expecting(PARAMETER))
211	{
212	expecting = Expecting(EQUALSIGN);
213	}
214	else
215	{
216	if (!strcmp(parmName, "Xxxxxxxxxx"))
217	{
218	expecting = Expecting(PARAMETER);
219	expecting \|= Expecting(SEMICOLON);
220	}
221	/* set Filename (plotParams.plotFileName) ***************************************************/
222	else if (!strcmp(parmName, "Filename"))
223	{
224	if (expecting == Expecting(EQUALSIGN))
225	{
226	expecting = Expecting(ALPHA);
227	readWord = YES;
228	}
229	else if (expecting == Expecting(ALPHA))
230	{
231	strcpy (plotParams.plotFileName, tkn);
232	MrBayesPrint ("%s Setting plot filename to %s\n", spacer, plotParams.plotFileName);
233	expecting = Expecting(PARAMETER) \| Expecting(SEMICOLON);
234	}
235	else
236	return (ERROR);
237	}
238	/* set Relburnin (chainParams.relativeBurnin) ********************************************************/
239	else if (!strcmp(parmName, "Relburnin"))
240	{
241	if (expecting == Expecting(EQUALSIGN))
242	expecting = Expecting(ALPHA);
243	else if (expecting == Expecting(ALPHA))
244	{
245	if (IsArgValid(tkn, tempStr) == NO_ERROR)
246	{
247	if (!strcmp(tempStr, "Yes"))
248	chainParams.relativeBurnin = YES;
249	else
250	chainParams.relativeBurnin = NO;
251	}
252	else
253	{
254	MrBayesPrint ("%s Invalid argument for Relburnin\n", spacer);
255	return (ERROR);
256	}
257	if (chainParams.relativeBurnin == YES)
258	MrBayesPrint ("%s Using relative burnin (a fraction of samples discarded).\n", spacer);
259	else
260	MrBayesPrint ("%s Using absolute burnin (a fixed number of samples discarded).\n", spacer);
261	expecting = Expecting(PARAMETER) \| Expecting(SEMICOLON);
262	}
263	else
264	{
265	return (ERROR);
266	}
267	}
268	/* set Burnin (chainParams.chainBurnIn) *******************************************************/
269	else if (!strcmp(parmName, "Burnin"))
270	{
271	if (expecting == Expecting(EQUALSIGN))
272	expecting = Expecting(NUMBER);
273	else if (expecting == Expecting(NUMBER))
274	{
275	sscanf (tkn, "%d", &tempI);
276	chainParams.chainBurnIn = tempI;
277	MrBayesPrint ("%s Setting burnin to %d\n", spacer, chainParams.chainBurnIn);
278	expecting = Expecting(PARAMETER) \| Expecting(SEMICOLON);
279	}
280	else
281	return (ERROR);
282	}
283	/* set Burninfrac (chainParams.burninFraction) ************************************************************/
284	else if (!strcmp(parmName, "Burninfrac"))
285	{
286	if (expecting == Expecting(EQUALSIGN))
287	expecting = Expecting(NUMBER);
288	else if (expecting == Expecting(NUMBER))
289	{
290	sscanf (tkn, "%lf", &tempD);
291	if (tempD < 0.01)
292	{
293	MrBayesPrint ("%s Burnin fraction too low (< 0.01)\n", spacer);
294	return (ERROR);
295	}
296	if (tempD > 0.50)
297	{
298	MrBayesPrint ("%s Burnin fraction too high (> 0.50)\n", spacer);
299	return (ERROR);
300	}
301	chainParams.burninFraction = tempD;
302	MrBayesPrint ("%s Setting burnin fraction to %.2f\n", spacer, chainParams.burninFraction);
303	expecting = Expecting(PARAMETER) \| Expecting(SEMICOLON);
304	}
305	else
306	{
307	return (ERROR);
308	}
309	}
310	/* set Parameter (plotParams.parameter) *******************************************************/
311	else if (!strcmp(parmName, "Parameter"))
312	{
313	if (expecting == Expecting(EQUALSIGN))
314	{
315	expecting = Expecting(ALPHA);
316	readWord = YES;
317	}
318	else if (expecting == Expecting(ALPHA))
319	{
320	strcpy (plotParams.parameter, tkn);
321	expecting = Expecting(PARAMETER) \| Expecting(SEMICOLON);
322	}
323	else
324	return (ERROR);
325	}
326	/* set Parameter (plotParams.match) *******************************************************/
327	else if (!strcmp(parmName, "Match"))
328	{
329	if (expecting == Expecting(EQUALSIGN))
330	expecting = Expecting(ALPHA);
331	else if (expecting == Expecting(ALPHA))
332	{
333	if (IsArgValid(tkn, tempStr) == NO_ERROR)
334	strcpy (plotParams.match, tempStr);
335	else
336	return (ERROR);
337
338	MrBayesPrint ("%s Setting plot matching to %s\n", spacer, plotParams.match);
339	expecting = Expecting(PARAMETER) \| Expecting(SEMICOLON);
340	}
341	else
342	return (ERROR);
343	}
344
345	else
346	return (ERROR);
347	}
348
349	return (NO_ERROR);
350
351	}
352

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

Context Navigation

source: trunk/GDE/MrBAYES/mrbayes_3.2.1/plot.c

Download in other formats: