Context Navigation

sump.c

Visit:

Last change on this file was 10416, checked in by aboeckma, 12 years ago
Added mr bayes (no menu yet)
File size: 93.1 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 "sump.h"
48	#include "mbmath.h"
49	#include "mcmc.h"
50	#include "utils.h"
51	#if defined(__MWERKS__)
52	#include "SIOUX.h"
53	#endif
54
55	const char* const svnRevisionSumpC="$Rev: 513 $"; /* Revision keyword which is expended/updated by svn on each commit/update*/
56
57	/* local prototypes */
58	int CompareModelProbs (const void x, const void y);
59	int PrintModelStats (char fileName, char headerNames, int nHeaders, ParameterSample parameterSamples, int nRuns, int nSamples);
60	int PrintOverlayPlot (MrBFlt xVals, MrBFlt yVals, int nRows, int startingFrom, int nSamples);
61	int PrintParamStats (char fileName, char headerNames, int nHeaders, ParameterSample parameterSamples, int nRuns, int nSamples);
62	void PrintPlotHeader (void);
63
64
65
66	/* AllocateParameterSamples: Allocate space for parameter samples */
67	int AllocateParameterSamples (ParameterSample **parameterSamples, int numRuns, int numRows, int numColumns)
68	{
69	int i, j;
70
71	(parameterSamples) = (ParameterSample ) SafeCalloc (numColumns, sizeof(ParameterSample));
72	if (!(*parameterSamples))
73	return ERROR;
74	(parameterSamples)[0].values = (MrBFlt ) SafeCalloc (numColumns numRuns, sizeof (MrBFlt *));
75	if (!((*parameterSamples)[0].values))
76	{
77	FreeParameterSamples(*parameterSamples);
78	return ERROR;
79	}
80	(parameterSamples)[0].values[0] = (MrBFlt ) SafeCalloc (numColumns * numRuns * numRows, sizeof (MrBFlt));
81	for (i=1; i<numColumns; i++)
82	(parameterSamples)[i].values = (parameterSamples)[0].values + i*numRuns;
83	for (i=1; i<numRuns; i++)
84	(parameterSamples)[0].values[i] = (parameterSamples)[0].values[0] + i*numRows;
85	for (i=1; i<numColumns; i++)
86	{
87	for (j=0; j<numRuns; j++)
88	{
89	(parameterSamples)[i].values[j] = (parameterSamples)[0].values[0] + inumRunsnumRows + j*numRows;
90	}
91	}
92
93	return NO_ERROR;
94	}
95
96
97
98
99
100	/** Compare function (ModelProb) for qsort. Note reverse sort order (from larger to smaller probs) */
101	int CompareModelProbs (const void x, const void y) {
102
103	if ((((ModelProb )(x))).prob > (((ModelProb )(y))).prob)
104	return -1;
105	else if ((((ModelProb )(x))).prob < (((ModelProb )(y))).prob)
106	return 1;
107	else
108	return 0;
109	}
110
111
112
113
114
115	int DoSump (void)
116
117	{
118
119	int i, n, nHeaders=0, numRows, numColumns, numRuns, whichIsX, whichIsY,
120	unreliable, oneUnreliable, burnin, longestHeader, len;
121	MrBFlt mean, harm_mean;
122	char **headerNames=NULL, temp[120];
123	SumpFileInfo fileInfo, firstFileInfo;
124	ParameterSample *parameterSamples=NULL;
125	FILE *fpLstat=NULL;
126
127
128	# if defined (MPI_ENABLED)
129	if (proc_id != 0)
130	return NO_ERROR;
131	# endif
132
133
134	/* tell user we are ready to go */
135	if (sumpParams.numRuns == 1)
136	MrBayesPrint ("%s Summarizing parameters in file %s.p\n", spacer, sumpParams.sumpFileName);
137	else if (sumpParams.numRuns == 2)
138	MrBayesPrint ("%s Summarizing parameters in files %s.run1.p and %s.run2.p\n", spacer, sumpParams.sumpFileName, sumpParams.sumpFileName);
139	else /* if (sumpParams.numRuns > 2) */
140	{
141	MrBayesPrint ("%s Summarizing parameters in %d files (%s.run1.p,\n", spacer, sumpParams.numRuns, sumpParams.sumpFileName);
142	MrBayesPrint ("%s %s.run2.p, etc)\n", spacer, sumpParams.sumpFileName);
143	}
144	MrBayesPrint ("%s Writing summary statistics to file %s.pstat\n", spacer, sumpParams.sumpFileName);
145
146	if (chainParams.relativeBurnin == YES)
147	MrBayesPrint ("%s Using relative burnin ('relburnin=yes'), discarding the first %.0f %% of samples\n",
148	spacer, chainParams.burninFraction*100.0, chainParams.burninFraction);
149	else
150	MrBayesPrint ("%s Using absolute burnin ('relburnin=no'), discarding the first %d samples\n",
151	spacer, chainParams.chainBurnIn, chainParams.chainBurnIn);
152
153	/* Initialize to silence warning. */
154	firstFileInfo.numRows = 0;
155	firstFileInfo.numColumns = 0;
156
157	/* examine input file(s) */
158	for (i=0; i<sumpParams.numRuns; i++)
159	{
160	if (sumpParams.numRuns == 1)
161	sprintf (temp, "%s.p", sumpParams.sumpFileName);
162	else
163	sprintf (temp, "%s.run%d.p", sumpParams.sumpFileName, i+1);
164
165	if (ExamineSumpFile (temp, &fileInfo, &headerNames, &nHeaders) == ERROR)
166	goto errorExit;
167
168	if (i==0)
169	{
170	if (fileInfo.numRows == 0 \|\| fileInfo.numColumns == 0)
171	{
172	MrBayesPrint ("%s The number of rows or columns in file %d is equal to zero\n", spacer, temp);
173	goto errorExit;
174	}
175	firstFileInfo = fileInfo;
176	}
177	else
178	{
179	if (firstFileInfo.numRows != fileInfo.numRows \|\| firstFileInfo.numColumns != fileInfo.numColumns)
180	{
181	MrBayesPrint ("%s First file had %d rows and %d columns while file %s had %d rows and %d columns\n",
182	spacer, firstFileInfo.numRows, firstFileInfo.numColumns, temp, fileInfo.numRows, fileInfo.numColumns);
183	MrBayesPrint ("%s MrBayes expects the same number of rows and columns in all files\n", spacer);
184	goto errorExit;
185	}
186	}
187	}
188
189	numRows = fileInfo.numRows;
190	numColumns = fileInfo.numColumns;
191	numRuns = sumpParams.numRuns;
192
193	/* allocate space to hold parameter information */
194	if (AllocateParameterSamples (&parameterSamples, numRuns, numRows, numColumns) == ERROR)
195	return ERROR;
196
197	/* read samples */
198	for (i=0; i<sumpParams.numRuns; i++)
199	{
200	/* derive file name */
201	if (sumpParams.numRuns == 1)
202	sprintf (temp, "%s.p", sumpParams.sumpFileName);
203	else
204	sprintf (temp, "%s.run%d.p", sumpParams.sumpFileName, i+1);
205
206	/* read samples */
207	if (ReadParamSamples (temp, &fileInfo, parameterSamples, i) == ERROR)
208	goto errorExit;
209	}
210
211	/* get length of longest header */
212	longestHeader = 9; /* 9 is the length of the word "parameter" (for printing table) */
213	for (i=0; i<nHeaders; i++)
214	{
215	len = (int) strlen(headerNames[i]);
216	if (len > longestHeader)
217	longestHeader = len;
218	}
219
220
221	/* Print header */
222	PrintPlotHeader ();
223
224	/* Print trace plots */
225	if (FindHeader("Gen", headerNames, nHeaders, &whichIsX) == ERROR)
226	{
227	MrBayesPrint ("%s Could not find the 'Gen' column\n", spacer);
228	return ERROR;
229	}
230	if (FindHeader("LnL", headerNames, nHeaders, &whichIsY) == ERROR)
231	{
232	MrBayesPrint ("%s Could not find the 'LnL' column\n", spacer);
233	return ERROR;
234	}
235
236
237	if (sumpParams.numRuns > 1)
238	{
239	if (sumpParams.allRuns == YES)
240	{
241	for (i=0; i<sumpParams.numRuns; i++)
242	{
243	MrBayesPrint ("\n%s Samples from run %d:\n", spacer, i+1);
244	if (PrintPlot (parameterSamples[whichIsX].values[i], parameterSamples[whichIsY].values[i], numRows) == ERROR)
245	goto errorExit;
246	}
247	}
248	else
249	{
250	if (PrintOverlayPlot (parameterSamples[whichIsX].values, parameterSamples[whichIsY].values, numRuns, 0, numRows) == ERROR)
251	goto errorExit;
252	}
253	}
254	else
255	{
256	if (PrintPlot (parameterSamples[whichIsX].values[0], parameterSamples[whichIsY].values[0], numRows) == ERROR)
257	goto errorExit;
258	}
259
260	/* calculate arithmetic and harmonic means of likelihoods */
261
262	/* open output file */
263	strncpy (temp, sumpParams.sumpOutfile, 90);
264	strcat (temp, ".lstat");
265	fpLstat = OpenNewMBPrintFile (temp);
266	if (!fpLstat)
267	goto errorExit;
268
269	/* print unique identifier to the output file */
270	if (strlen(stamp) > 1)
271	fprintf (fpLstat, "[ID: %s]\n", stamp);
272
273	/* print header */
274	if (sumpParams.numRuns == 1)
275	MrBayesPrintf(fpLstat, "arithmetic_mean\tharmonic_mean\tvalues_discarded\n");
276	else
277	MrBayesPrintf(fpLstat, "run\tarithmetic_mean\tharmonic_mean\tvalues_discarded\n");
278
279	oneUnreliable = NO;
280	for (n=0; n<sumpParams.numRuns; n++)
281	{
282	unreliable = NO;
283	if (HarmonicArithmeticMeanOnLogs (parameterSamples[whichIsY].values[n], numRows, &mean, &harm_mean) == ERROR)
284	{
285	unreliable = YES;
286	oneUnreliable = YES;
287	}
288	if (sumpParams.numRuns == 1)
289	{
290	MrBayesPrint ("\n");
291	MrBayesPrint ("%s Estimated marginal likelihoods for run sampled in file \"%s.p\":\n", spacer, sumpParams.sumpFileName);
292	MrBayesPrint ("%s (Use the harmonic mean for Bayes factor comparisons of models)\n", spacer, sumpParams.sumpFileName);
293	MrBayesPrint ("%s (Values are saved to the file %s.lstat)\n\n", spacer, sumpParams.sumpOutfile);
294	MrBayesPrint ("%s Arithmetic mean Harmonic mean\n", spacer);
295	MrBayesPrint ("%s --------------------------------\n", spacer);
296	if (unreliable == NO)
297	MrBayesPrint ("%s %9.2lf %9.2lf\n", spacer, mean, harm_mean);
298	else
299	MrBayesPrint ("%s %9.2lf * %9.2lf *\n", spacer, mean, harm_mean);
300
301	/* print to file */
302	MrBayesPrintf(fpLstat, "%s\t", MbPrintNum(mean));
303	MrBayesPrintf(fpLstat, "%s\t", MbPrintNum(harm_mean));
304	if (unreliable == YES)
305	MrBayesPrintf(fpLstat, "yes\n");
306	else
307	MrBayesPrintf(fpLstat, "no\n");
308	}
309	else
310	{
311	if (n == 0)
312	{
313	MrBayesPrint ("\n");
314	MrBayesPrint ("%s Estimated marginal likelihoods for runs sampled in files\n", spacer);
315	if (sumpParams.numRuns > 2)
316	MrBayesPrint ("%s \"%s.run1.p\", \"%s.run2.p\", etc:\n", spacer, sumpParams.sumpFileName, sumpParams.sumpFileName);
317	else /* if (sumpParams.numRuns == 2) */
318	MrBayesPrint ("%s \"%s.run1.p\" and \"%s.run2.p\":\n", spacer, sumpParams.sumpFileName, sumpParams.sumpFileName);
319	MrBayesPrint ("%s (Use the harmonic mean for Bayes factor comparisons of models)\n\n", spacer, sumpParams.sumpFileName);
320	MrBayesPrint ("%s (Values are saved to the file %s.lstat)\n\n", spacer, sumpParams.sumpOutfile);
321	MrBayesPrint ("%s Run Arithmetic mean Harmonic mean\n", spacer);
322	MrBayesPrint ("%s --------------------------------------\n", spacer);
323	}
324	if (unreliable == NO)
325	MrBayesPrint ("%s %3d %9.2lf %9.2lf\n", spacer, n+1, mean, harm_mean);
326	else
327	MrBayesPrint ("%s %3d %9.2lf * %9.2lf *\n", spacer, n+1, mean, harm_mean);
328
329	/* print to file */
330	MrBayesPrintf(fpLstat, "%d\t", n+1);
331	MrBayesPrintf(fpLstat, "%s\t", MbPrintNum(mean));
332	MrBayesPrintf(fpLstat, "%s\t", MbPrintNum(harm_mean));
333	if (unreliable == YES)
334	MrBayesPrintf(fpLstat, "yes\n");
335	else
336	MrBayesPrintf(fpLstat, "no\n");
337	}
338	} /* next run */
339	if (sumpParams.numRuns == 1)
340	{
341	MrBayesPrint ("%s --------------------------------\n", spacer);
342	}
343	else
344	{
345	if (HarmonicArithmeticMeanOnLogs (parameterSamples[whichIsY].values[0], sumpParams.numRuns*numRows, &mean, &harm_mean) == ERROR)
346	{
347	unreliable = YES;
348	oneUnreliable = YES;
349	}
350	else
351	unreliable = NO;
352	MrBayesPrint ("%s --------------------------------------\n", spacer);
353	if (unreliable == YES)
354	MrBayesPrint ("%s TOTAL %9.2lf * %9.2lf *\n", spacer, mean, harm_mean);
355	else
356	MrBayesPrint ("%s TOTAL %9.2lf %9.2lf\n", spacer, mean, harm_mean);
357	MrBayesPrint ("%s --------------------------------------\n", spacer);
358
359	/* print total to file */
360	MrBayesPrintf(fpLstat, "all\t");
361	MrBayesPrintf(fpLstat, "%s\t", MbPrintNum(mean));
362	MrBayesPrintf(fpLstat, "%s\t", MbPrintNum(harm_mean));
363	if (unreliable == YES)
364	MrBayesPrintf(fpLstat, "yes\n");
365	else
366	MrBayesPrintf(fpLstat, "no\n");
367	}
368	if (oneUnreliable == YES)
369	{
370	MrBayesPrint ("%s * These estimates may be unreliable because \n", spacer);
371	MrBayesPrint ("%s some extreme values were excluded\n\n", spacer);
372	}
373	else
374	{
375	MrBayesPrint ("\n");
376	}
377	SafeFclose(&fpLstat);
378
379	/* Calculate burnin */
380	burnin = fileInfo.firstParamLine - fileInfo.headerLine;
381
382	/* Print parameter information to screen and to file. */
383	if (sumpParams.numRuns > 1 && sumpParams.allRuns == YES)
384	{
385	for (i=0; i<sumpParams.numRuns; i++)
386	{
387	/* print table header */
388	MrBayesPrint ("\n");
389	MrBayesPrint ("%s Model parameter summaries for run sampled in file \"%s.run%d.p\":\n", spacer, sumpParams.sumpFileName, i+1);
390	MrBayesPrint ("%s (Based on %d samples out of a total of %d samples from this analysis)\n\n", spacer, numRows, numRows + burnin);
391	if (PrintParamStats (sumpParams.sumpOutfile, headerNames, nHeaders, parameterSamples, numRuns, numRows) == ERROR)
392	goto errorExit;
393	if (PrintModelStats (sumpParams.sumpOutfile, headerNames, nHeaders, parameterSamples, numRuns, numRows) == ERROR)
394	goto errorExit;
395	}
396	}
397
398	MrBayesPrint ("\n");
399	if (sumpParams.numRuns == 1)
400	MrBayesPrint ("%s Model parameter summaries for run sampled in file \"%s\":\n", spacer, sumpParams.sumpFileName);
401	else if (sumpParams.numRuns == 2)
402	{
403	MrBayesPrint ("%s Model parameter summaries over the runs sampled in files\n", spacer);
404	MrBayesPrint ("%s \"%s.run1.p\" and \"%s.run2.p\":\n", spacer, sumpParams.sumpFileName, sumpParams.sumpFileName);
405	}
406	else
407	{
408	MrBayesPrint ("%s Model parameter summaries over all %d runs sampled in files\n", spacer, sumpParams.numRuns);
409	MrBayesPrint ("%s \"%s.run1.p\", \"%s.run2.p\" etc:\n", spacer, sumpParams.sumpFileName, sumpParams.sumpFileName);
410	}
411
412	if (sumpParams.numRuns == 1)
413	{
414	MrBayesPrint ("%s Based on a total of %d samples out of a total of %d samples\n", spacer, numRows, numRows + burnin);
415	MrBayesPrint ("%s from this analysis.\n", spacer);
416	}
417	else
418	{
419	MrBayesPrint ("%s Summaries are based on a total of %d samples from %d runs.\n", spacer, sumpParams.numRuns*numRows, sumpParams.numRuns);
420	MrBayesPrint ("%s Each run produced %d samples of which %d samples were included.\n", spacer, numRows + burnin, numRows);
421	}
422	MrBayesPrint ("%s Parameter summaries saved to file \"%s.pstat\".\n", spacer, sumpParams.sumpOutfile);
423
424	if (PrintParamStats (sumpParams.sumpOutfile, headerNames, nHeaders, parameterSamples, numRuns, numRows) == ERROR)
425	goto errorExit;
426	if (PrintModelStats (sumpParams.sumpOutfile, headerNames, nHeaders, parameterSamples, numRuns, numRows) == ERROR)
427	goto errorExit;
428
429	/* free memory */
430	FreeParameterSamples(parameterSamples);
431	for (i=0; i<nHeaders; i++)
432	free (headerNames[i]);
433	free (headerNames);
434
435	expecting = Expecting(COMMAND);
436	strcpy (spacer, "");
437
438	return (NO_ERROR);
439
440	errorExit:
441
442	/* free memory */
443	FreeParameterSamples (parameterSamples);
444	for (i=0; i<nHeaders; i++)
445	free (headerNames[i]);
446	free (headerNames);
447
448	if (fpLstat)
449	SafeFclose (&fpLstat);
450
451	expecting = Expecting(COMMAND);
452	strcpy (spacer, "");
453
454	return (ERROR);
455	}
456
457
458
459
460	int DoSumSs (void)
461	{
462
463	int i, nHeaders=0, numRows, numColumns, numRuns, whichIsX, whichIsY,
464	longestHeader, len;
465	char **headerNames=NULL, temp[120];
466	SumpFileInfo fileInfo, firstFileInfo;
467	ParameterSample *parameterSamples=NULL;
468	int stepIndexSS,numSamplesInStepSS, stepBeginSS, stepBurnin;
469	MrBFlt lnlp, nextSteplnlp, *firstlnlp;
470	MrBFlt *marginalLnLSS=NULL,stepScalerSS,stepAcumulatorSS, stepLengthSS, tmpMfl;
471	int beginPrint, countPrint;
472	float tmpf;
473	MrBFlt plotArrayY=NULL,plotArrayX=NULL;
474	int j,k,count,result;
475	MrBFlt sum;
476	int firstPass = YES;
477
478	# if defined (MPI_ENABLED)
479	if (proc_id != 0)
480	return NO_ERROR;
481	# endif
482
483	chainParams.isSS=YES;
484
485	/* tell user we are ready to go */
486	if (sumssParams.numRuns == 1)
487	MrBayesPrint ("%s Summarizing parameters in file %s.p\n", spacer, sumpParams.sumpFileName);
488	else if (sumssParams.numRuns == 2)
489	MrBayesPrint ("%s Summarizing parameters in files %s.run1.p and %s.run2.p\n", spacer, sumpParams.sumpFileName, sumpParams.sumpFileName);
490	else /* if (sumssParams.numRuns > 2) */
491	{
492	MrBayesPrint ("%s Summarizing parameters in %d files (%s.run1.p,\n", spacer, sumssParams.numRuns, sumpParams.sumpFileName);
493	MrBayesPrint ("%s %s.run2.p, etc)\n", spacer, sumpParams.sumpFileName);
494	}
495	//MrBayesPrint ("%s Writing summary statistics to file %s.pstat\n", spacer, sumpParams.sumpFileName);
496
497	if (chainParams.relativeBurnin == YES)
498	MrBayesPrint ("%s Using relative burnin ('relburnin=yes'), discarding the first %.0f %% of samples within each step.\n",
499	spacer, chainParams.burninFraction*100.0, chainParams.burninFraction);
500	else
501	MrBayesPrint ("%s Using absolute burnin ('relburnin=no'), discarding the first %d samples within each step.\n",
502	spacer, chainParams.chainBurnIn, chainParams.chainBurnIn);
503
504	/* Initialize to silence warning. */
505	firstFileInfo.numRows = 0;
506	firstFileInfo.numColumns = 0;
507
508	/* examine input file(s) */
509	for (i=0; i<sumssParams.numRuns; i++)
510	{
511	if (sumssParams.numRuns == 1)
512	sprintf (temp, "%s.p", sumpParams.sumpFileName);
513	else
514	sprintf (temp, "%s.run%d.p", sumpParams.sumpFileName, i+1);
515
516	if (ExamineSumpFile (temp, &fileInfo, &headerNames, &nHeaders) == ERROR)
517	goto errorExit;
518
519	if (i==0)
520	{
521	if (fileInfo.numRows == 0 \|\| fileInfo.numColumns == 0)
522	{
523	MrBayesPrint ("%s The number of rows or columns in file %d is equal to zero\n", spacer, temp);
524	goto errorExit;
525	}
526	firstFileInfo = fileInfo;
527	}
528	else
529	{
530	if (firstFileInfo.numRows != fileInfo.numRows \|\| firstFileInfo.numColumns != fileInfo.numColumns)
531	{
532	MrBayesPrint ("%s First file had %d rows and %d columns while file %s had %d rows and %d columns\n",
533	spacer, firstFileInfo.numRows, firstFileInfo.numColumns, temp, fileInfo.numRows, fileInfo.numColumns);
534	MrBayesPrint ("%s MrBayes expects the same number of rows and columns in all files\n", spacer);
535	goto errorExit;
536	}
537	}
538	}
539
540	numRows = fileInfo.numRows;
541	numColumns = fileInfo.numColumns;
542	numRuns = sumssParams.numRuns;
543
544	/* allocate space to hold parameter information */
545	if (AllocateParameterSamples (&parameterSamples, numRuns, numRows, numColumns) == ERROR)
546	goto errorExit;
547
548	/* read samples */
549	for (i=0; i<sumssParams.numRuns; i++)
550	{
551	/* derive file name */
552	if (sumssParams.numRuns == 1)
553	sprintf (temp, "%s.p", sumpParams.sumpFileName);
554	else
555	sprintf (temp, "%s.run%d.p", sumpParams.sumpFileName, i+1);
556
557	/* read samples */
558	if (ReadParamSamples (temp, &fileInfo, parameterSamples, i) == ERROR)
559	goto errorExit;
560	}
561
562	/* get length of longest header */
563	longestHeader = 9; /* 9 is the length of the word "parameter" (for printing table) */
564	for (i=0; i<nHeaders; i++)
565	{
566	len = (int) strlen(headerNames[i]);
567	if (len > longestHeader)
568	longestHeader = len;
569	}
570
571
572	/* Print trace plots */
573	if (FindHeader("Gen", headerNames, nHeaders, &whichIsX) == ERROR)
574	{
575	MrBayesPrint ("%s Could not find the 'Gen' column\n", spacer);
576	goto errorExit;
577	}
578	if (FindHeader("LnL", headerNames, nHeaders, &whichIsY) == ERROR)
579	{
580	MrBayesPrint ("%s Could not find the 'LnL' column\n", spacer);
581	goto errorExit;
582	}
583
584
585
586	if(chainParams.burninSS > 0)
587	{
588	stepBeginSS = chainParams.burninSS + 1;
589	}
590	else
591	{
592	numSamplesInStepSS = (numRows-1)/(chainParams.numStepsSS-chainParams.burninSS);
593	stepBeginSS = numSamplesInStepSS + 1;
594	}
595
596	numSamplesInStepSS = (numRows - stepBeginSS)/chainParams.numStepsSS;
597	if( (numRows - stepBeginSS)%chainParams.numStepsSS!=0 )
598	{
599	MrBayesPrint ("%s Error: Number of samples could not be evenly devided among steps (%d samples among %d steps). \n", spacer,(numRows - stepBeginSS),chainParams.numStepsSS);
600	goto errorExit;
601	}
602
603
604	if( chainParams.relativeBurnin == YES )
605	{
606	stepBurnin = (int)(numSamplesInStepSS*chainParams.burninFraction);
607	}
608	else
609	{
610	stepBurnin = chainParams.chainBurnIn;
611	if(stepBurnin >= numSamplesInStepSS )
612	{
613	MrBayesPrint ("%s Error: Burnin in each step(%d) is longer then the step itself(%d). \n", spacer,stepBurnin, numSamplesInStepSS );
614	goto errorExit;
615	}
616	}
617
618	marginalLnLSS = (MrBFlt *) SafeCalloc (sumssParams.numRuns, sizeof(MrBFlt));
619	/Preparing and printing joined plot./
620	plotArrayY = (MrBFlt *) SafeCalloc (sumssParams.numRuns+1, sizeof(MrBFlt));
621	for(i=0; i<sumssParams.numRuns+1; i++)
622	plotArrayY[i] = (MrBFlt *) SafeCalloc (numSamplesInStepSS, sizeof(MrBFlt));
623
624	plotArrayX = (MrBFlt *) SafeCalloc (sumssParams.numRuns, sizeof(MrBFlt));
625	for(i=0; i<sumssParams.numRuns; i++)
626	{
627	plotArrayX[i] = (MrBFlt *) SafeCalloc (numSamplesInStepSS, sizeof(MrBFlt));
628	for(j=0; j<numSamplesInStepSS; j++)
629	plotArrayX[i][j]=j+1;
630	}
631
632	MrBayesPrint ("%s In total %d sampls are red from .p files.\n", spacer, numRows );
633	MrBayesPrint ("\n");
634	MrBayesPrint ("%s Marginal likelihood (in natural log units) is estimated using stepping-stone sampling\n", spacer );
635	MrBayesPrint ("%s based on %d steps with %d samples within each step. \n", spacer, chainParams.numStepsSS, numSamplesInStepSS );
636	MrBayesPrint ("%s First %d samples (including generation 0) are discarded as initial burn-in.\n", spacer, stepBeginSS);
637	if(chainParams.startFromPriorSS==YES)
638	MrBayesPrint ("%s Sampling is assumed have being done from prior to posterior.\n", spacer);
639	else
640	{
641	MrBayesPrint ("%s Sampling is assumed have being done from posterior to prior.\n", spacer);
642	}
643
644	sumssTable:
645
646	MrBayesPrint ("\n\n%s Step contribution table.\n\n",spacer);
647	MrBayesPrint (" Columns in the table: \n");
648	MrBayesPrint (" Step -- Index of the step \n");
649	MrBayesPrint (" runX -- Contribution to the marginal log likelihood of run X, i.e. marginal \n");
650	MrBayesPrint (" log likelihood for run X is the sum across all steps in column runX.\n\n");
651
652	if( firstPass == YES && chainParams.relativeBurnin == YES )
653	MrBayesPrint ("%s The table entrances are based on samples excluding burn-in %d samples (%d%%) \n", spacer, stepBurnin,(int)(100*chainParams.burninFraction) );
654	else
655	MrBayesPrint ("%s The table entrances are based on samples excluding burn-in %d samples \n", spacer, stepBurnin);
656	MrBayesPrint ("%s discarded at the begining of each step. \n\n", spacer);
657
658	//MrBayesPrint ("%s Run Marginal likelihood (ln)\n",spacer);
659	//MrBayesPrint ("%s ------------------------------\n",spacer);
660	MrBayesPrint (" Step");
661	for (j=0; j<sumssParams.numRuns ; j++)
662	{
663	if(j<9)
664	MrBayesPrint (" ");
665	MrBayesPrint (" run%d", j+1);
666	}
667	MrBayesPrint ("\n");
668	for(i=0; i<sumssParams.numRuns; i++)
669	{
670	marginalLnLSS[i] = 0.0;
671	}
672	for(stepIndexSS = chainParams.numStepsSS-1; stepIndexSS>=0; stepIndexSS--)
673	{
674	if(chainParams.startFromPriorSS==YES)
675	{
676	stepLengthSS = BetaQuantile( chainParams.alphaSS, 1.0, (MrBFlt)(chainParams.numStepsSS-stepIndexSS)/(MrBFlt)chainParams.numStepsSS)-BetaQuantile( chainParams.alphaSS, 1.0, (MrBFlt)(chainParams.numStepsSS-1-stepIndexSS)/(MrBFlt)chainParams.numStepsSS);
677	}
678	else
679	{
680	stepLengthSS = BetaQuantile ( chainParams.alphaSS, 1.0, (MrBFlt)(stepIndexSS+1)/(MrBFlt)chainParams.numStepsSS) - BetaQuantile ( chainParams.alphaSS, 1.0, (MrBFlt)stepIndexSS/(MrBFlt)chainParams.numStepsSS);
681	}
682	MrBayesPrint (" %3d ", chainParams.numStepsSS-stepIndexSS);
683	for(i=0; i<sumssParams.numRuns; i++)
684	{
685	lnlp = parameterSamples[whichIsY].values[i] + stepBeginSS + (chainParams.numStepsSS-stepIndexSS-1)*numSamplesInStepSS;
686	nextSteplnlp = lnlp+numSamplesInStepSS;
687	lnlp+= stepBurnin;
688	stepAcumulatorSS = 0.0;
689	stepScalerSS = lnlpstepLengthSS;
690	while( lnlp<nextSteplnlp )
691	{
692	if( lnlpstepLengthSS > stepScalerSS + 200.0 )
693	{
694	// adjust scaler;
695	stepAcumulatorSS /= exp( lnlpstepLengthSS - 100.0 - stepScalerSS );
696	stepScalerSS= lnlpstepLengthSS - 100.0;
697	}
698	stepAcumulatorSS += exp( lnlpstepLengthSS - stepScalerSS );
699	lnlp++;
700	}
701	tmpMfl = (log( stepAcumulatorSS/(numSamplesInStepSS-stepBurnin) ) + stepScalerSS);
702	MrBayesPrint (" %10.3lf", tmpMfl);
703	marginalLnLSS[i] += tmpMfl;
704	}
705	MrBayesPrint ("\n");
706	//MrBayesPrint ("%s %3d %9.2f \n", spacer, i+1, marginalLnLSS );
707	}
708	MrBayesPrint (" ");
709	for (j=0; j<sumssParams.numRuns ; j++)
710	{
711	if(j<9)
712	MrBayesPrint ("-");
713	MrBayesPrint ("----------");
714	}
715	MrBayesPrint ("\n");
716	MrBayesPrint (" Sum: ");
717	for (j=0; j<sumssParams.numRuns ; j++)
718	MrBayesPrint (" %10.3lf", marginalLnLSS[j]);
719
720	MrBayesPrint ("\n");
721	/*
722	if (sumssParams.numRuns > 1)
723	{
724	MrBayesPrint ("%s Below are rough plots of the generations (x-axis) during burn in \n", spacer);
725	MrBayesPrint ("%s phase versus the log probability of observing the data (y-axis). \n", spacer);
726	MrBayesPrint ("%s You can use these graphs to determine if the burn in for your SS \n", spacer);
727	MrBayesPrint ("%s analysis was sufficiant. The log probability suppose to plateau \n", spacer);
728	MrBayesPrint ("%s indicating that you may be at stationarity by the time you finish \n", spacer);
729	MrBayesPrint ("%s burn in phase. This burn in, unlike burn in within each step, is \n", spacer);
730	MrBayesPrint ("%s fixed and can not be changed. \n", spacer);
731	}
732	else
733	{
734	MrBayesPrint ("%s Below is a rough plot of the generations (x-axis) during burn in \n", spacer);
735	MrBayesPrint ("%s phase versus the log probability of observing the data (y-axis). \n", spacer);
736	MrBayesPrint ("%s You can use these graph to determine if the burn in for your SS \n", spacer);
737	MrBayesPrint ("%s analysis was sufficiant. The log probability suppose to plateau \n", spacer);
738	MrBayesPrint ("%s indicating that you may be at stationarity by the time you finish \n", spacer);
739	MrBayesPrint ("%s burn in phase. This burn in, unlike burn in within each step, is \n", spacer);
740	MrBayesPrint ("%s fixed and can not be changed. \n", spacer);
741	}
742	*/
743
744	if( firstPass == NO )
745	goto sumssExitOptions;
746
747	sumssStepPlot:
748
749	MrBayesPrint ("\n\n%s Step plot(s).\n",spacer);
750	while(1)
751	{
752	MrBayesPrint ("\n");
753	if( sumssParams.stepToPlot == 0 )
754	{
755	beginPrint=(int)(sumssParams.discardFraction*stepBeginSS);
756	countPrint=stepBeginSS-beginPrint;
757	MrBayesPrint ("%s Ploting step 0, i.e initial burn-in phase consisting of %d samples.\n", spacer,stepBeginSS);
758	MrBayesPrint ("%s According to 'Discardfrac=%.2f', first %d samples are not ploted.\n", spacer,sumssParams.discardFraction,beginPrint);
759	}
760	else
761	{
762	if( sumssParams.stepToPlot > chainParams.numStepsSS )
763	{
764	MrBayesPrint ("%s Chosen index of step to print %d is out of range of step indices[0,...,%d].\n", spacer,sumssParams.stepToPlot,chainParams.numStepsSS);
765	goto errorExit;
766	}
767	beginPrint=stepBeginSS+(sumssParams.stepToPlot-1)numSamplesInStepSS + (int)(sumssParams.discardFractionnumSamplesInStepSS);
768	countPrint=numSamplesInStepSS-(int)(sumssParams.discardFraction*numSamplesInStepSS);
769	MrBayesPrint ("%s Ploting step %d consisting of %d samples.\n", spacer,sumssParams.stepToPlot,numSamplesInStepSS);
770	MrBayesPrint ("%s According to 'Discardfrac=%.2f', first %d samples are not ploted.\n", spacer,sumssParams.discardFraction,(int)(sumssParams.discardFraction*numSamplesInStepSS));
771	}
772
773
774	if (sumssParams.numRuns > 1)
775	{
776	if (sumpParams.allRuns == YES)
777	{
778	for (i=0; i<sumssParams.numRuns; i++)
779	{
780	MrBayesPrint ("\n%s Samples from run %d:\n", spacer, i+1);
781	if (PrintPlot (parameterSamples[whichIsX].values[i]+beginPrint, parameterSamples[whichIsY].values[i]+beginPrint, countPrint) == ERROR)
782	goto errorExit;
783	}
784	}
785	else
786	{
787	if (PrintOverlayPlot (parameterSamples[whichIsX].values, parameterSamples[whichIsY].values, numRuns, beginPrint, countPrint) == ERROR)
788	goto errorExit;
789	}
790	}
791	else
792	{
793	if (PrintPlot (parameterSamples[whichIsX].values[0]+beginPrint, parameterSamples[whichIsY].values[0]+beginPrint, countPrint) == ERROR)
794	goto errorExit;
795	}
796
797	if( sumssParams.askForMorePlots == NO \|\| firstPass == YES )
798	break;
799
800	MrBayesPrint (" You can choose to print new step plots for different steps or discard fractions.\n");
801	MrBayesPrint (" Allowed range of 'Steptoplot' are from 0 to %d.\n", chainParams.numStepsSS);
802	MrBayesPrint (" If the next entered value is negative, 'sumss' will stop printing step plots.\n");
803	MrBayesPrint (" If the next entered value is positive, but out of range, you will be offered\n");
804	MrBayesPrint (" to change paramiter 'Discardfrac' of 'sumss'.\n");
805	MrBayesPrint (" Enter new step number 'Steptoplot':");
806	result=scanf("%d",&j);
807	if(j < 0 )
808	break;
809	if(j > chainParams.numStepsSS)
810	{
811	do
812	{
813	MrBayesPrint (" Enter new value for 'Discardfrac', should be in range 0.0 to 1.0:");
814	result=scanf("%f",&tmpf);
815	sumssParams.discardFraction = (MrBFlt)tmpf;
816	}
817	while(sumssParams.discardFraction < 0.0 \|\| sumssParams.discardFraction > 1.0);
818	}
819	else
820	sumssParams.stepToPlot=j;
821	}
822
823	if( firstPass == NO )
824	goto sumssExitOptions;
825
826	sumssJoinedPlot:
827
828	MrBayesPrint ("\n\n%s Joined plot(s).\n",spacer);
829	while(1)
830	{
831	MrBayesPrint ("\n");
832	MrBayesPrint ("%s Joined plot of %d samples of all steps together. 'smoothing' is set to:%d\n", spacer,numSamplesInStepSS,sumssParams.smoothing);
833	MrBayesPrint ("%s According to step burn-in, first %d samples are not ploted.\n", spacer,stepBurnin);
834
835	for(i=0; i<sumssParams.numRuns; i++)
836	{
837	for(j=stepBurnin;j<numSamplesInStepSS;j++)
838	plotArrayY[sumssParams.numRuns][j]=0.0;
839	lnlp= parameterSamples[whichIsY].values[i] + stepBeginSS;
840	nextSteplnlp=lnlp;
841	for(stepIndexSS = chainParams.numStepsSS-1; stepIndexSS>0; stepIndexSS--)
842	{
843	firstlnlp=plotArrayY[sumssParams.numRuns] + stepBurnin;
844	lnlp+=stepBurnin;
845	nextSteplnlp += numSamplesInStepSS;
846	while( lnlp<nextSteplnlp )
847	{
848	firstlnlp+=lnlp;
849	firstlnlp++;
850	lnlp++;
851	}
852	}
853	for(j=stepBurnin;j<numSamplesInStepSS;j++)
854	{
855	sum=0.0;
856	count=0;
857	for(k=j-sumssParams.smoothing;k<=j+sumssParams.smoothing;k++)
858	{
859	if(k>=stepBurnin && k<numSamplesInStepSS)
860	{
861	sum += plotArrayY[sumssParams.numRuns][k];
862	count++;
863	}
864	}
865	plotArrayY[i][j] = sum/count;
866	/*
867	if( max < plotArrayY[i][j])
868	max=plotArrayY[i][j];
869	*/
870	}
871	/* for(j=stepBurnin;j<numSamplesInStepSS;j++)
872	{
873	plotArrayY[i][j] /= max;
874	}*/
875	}
876
877	beginPrint=stepBurnin;
878	countPrint=numSamplesInStepSS-stepBurnin;
879
880	if (sumssParams.numRuns > 1)
881	{
882	if (sumpParams.allRuns == YES)
883	{
884	for (i=0; i<sumssParams.numRuns; i++)
885	{
886	MrBayesPrint ("\n%s Samples from run %d:\n", spacer, i+1);
887	if (PrintPlot (plotArrayX[i]+beginPrint, plotArrayY[i]+beginPrint, countPrint) == ERROR)
888	goto errorExit;
889	}
890	}
891	else
892	{
893	if (PrintOverlayPlot (plotArrayX, plotArrayY, numRuns, beginPrint, countPrint) == ERROR)
894	goto errorExit;
895	}
896	}
897	else
898	{
899	if (PrintPlot (plotArrayX[0]+beginPrint, plotArrayY[0]+beginPrint, countPrint) == ERROR)
900	goto errorExit;
901	}
902
903	if( sumssParams.askForMorePlots == NO \|\| firstPass == YES )
904	break;
905
906	MrBayesPrint (" You can choose to print new joined plots with different step burn-in or smoothing.\n");
907	MrBayesPrint (" Allowed range of step burn-in values are from 0 to %d.\n", numSamplesInStepSS-1);
908	MrBayesPrint (" If the next entered value is negative, 'sumss' will stop printing joined plots.\n");
909	MrBayesPrint (" If the next entered value is positive, but out of range, you will be offered\n");
910	MrBayesPrint (" to change 'Smoothimg'.\n");
911	MrBayesPrint (" Enter new step burn-in:");
912	result=scanf("%d",&j);
913	if(j < 0 )
914	break;
915	if(j >= numSamplesInStepSS)
916	{
917	MrBayesPrint (" Enter new value for 'Smoothing':");
918	result=scanf("%d",&j);
919	sumssParams.smoothing = abs(j);
920	}
921	else
922	stepBurnin=j;
923	}
924
925	firstPass = NO;
926	sumssExitOptions:
927	if(sumssParams.askForMorePlots == YES )
928	{
929	MrBayesPrint ("\n");
930	MrBayesPrint (" Sumss is interactive, because of paramiter 'Askmore=YES' setting. \n");
931	MrBayesPrint (" What would you like to do next?\n");
932	MrBayesPrint (" 1) Print updated table according to new step burn-in.\n");
933	MrBayesPrint (" 2) Print Step plot(s).\n");
934	MrBayesPrint (" 3) Print Joined plot(s).\n");
935	MrBayesPrint (" 4) Exit 'sumss'.\n");
936	MrBayesPrint (" Enter a number that corresponds to one of the options:");
937	do
938	{
939	result=scanf("%d",&j);
940	}while(j<1 \|\| j>4);
941
942	if(j == 1)
943	{
944	MrBayesPrint (" Allowed range of step burn-in values are from 0 to %d\n", numSamplesInStepSS-1);
945	MrBayesPrint (" Current step burn-in value is:%d\n", stepBurnin);
946	MrBayesPrint (" Enter new step burn-in:");
947	do
948	{
949	result=scanf("%d",&stepBurnin);
950	}
951	while(stepBurnin < 0 \|\| stepBurnin > numSamplesInStepSS-1);
952	MrBayesPrint ("\n");
953	goto sumssTable;
954	}
955	else if(j == 2)
956	{
957	goto sumssStepPlot;
958	}
959	else if(j == 3)
960	goto sumssJoinedPlot;
961
962	}
963
964	/* free memory */
965	FreeParameterSamples(parameterSamples);
966	for (i=0; i<nHeaders; i++)
967	free (headerNames[i]);
968	free (headerNames);
969
970	expecting = Expecting(COMMAND);
971	strcpy (spacer, "");
972	chainParams.isSS=NO;
973	for(i=0; i<sumssParams.numRuns+1; i++)
974	free(plotArrayY[i]);
975	free(plotArrayY);
976	for(i=0; i<sumssParams.numRuns; i++)
977	free(plotArrayX[i]);
978	free(plotArrayX);
979	free(marginalLnLSS);
980
981	return (NO_ERROR);
982
983	errorExit:
984
985	/* free memory */
986	FreeParameterSamples (parameterSamples);
987	if( headerNames!=NULL )
988	for (i=0; i<nHeaders; i++)
989	free (headerNames[i]);
990	free (headerNames);
991
992	expecting = Expecting(COMMAND);
993	strcpy (spacer, "");
994	chainParams.isSS=NO;
995	if( plotArrayY!=NULL )
996	for(i=0; i<sumssParams.numRuns+1; i++)
997	free(plotArrayY[i]);
998	free(plotArrayY);
999	if( plotArrayX!=NULL )
1000	for(i=0; i<sumssParams.numRuns; i++)
1001	free(plotArrayX[i]);
1002	free(plotArrayX);
1003	free(marginalLnLSS);
1004
1005	return (ERROR);
1006	}
1007
1008
1009
1010
1011
1012	int DoSumpParm (char parmName, char tkn)
1013
1014	{
1015
1016	int tempI;
1017	MrBFlt tempD;
1018	char tempStr[100];
1019
1020	if (expecting == Expecting(PARAMETER))
1021	{
1022	expecting = Expecting(EQUALSIGN);
1023	}
1024	else
1025	{
1026	if (!strcmp(parmName, "Xxxxxxxxxx"))
1027	{
1028	expecting = Expecting(PARAMETER);
1029	expecting \|= Expecting(SEMICOLON);
1030	}
1031	/* set Filename (sumpParams.sumpFileName) ***************************************************/
1032	else if (!strcmp(parmName, "Filename"))
1033	{
1034	if (expecting == Expecting(EQUALSIGN))
1035	{
1036	expecting = Expecting(ALPHA);
1037	readWord = YES;
1038	}
1039	else if (expecting == Expecting(ALPHA))
1040	{
1041	if(strlen(tkn)>99 && (strchr(tkn,' ')-tkn) > 99 )
1042	{
1043	MrBayesPrint ("%s Maximum allowed length of file name is 99 characters. The given name:\n", spacer);
1044	MrBayesPrint ("%s '%s'\n", spacer,tkn);
1045	return (ERROR);
1046	}
1047	sscanf (tkn, "%s", tempStr);
1048	strcpy (sumpParams.sumpFileName, tempStr);
1049	strcpy (sumpParams.sumpOutfile, tempStr);
1050	MrBayesPrint ("%s Setting sump filename and output file name to %s\n", spacer, sumpParams.sumpFileName);
1051	expecting = Expecting(PARAMETER) \| Expecting(SEMICOLON);
1052	}
1053	else
1054	return (ERROR);
1055	}
1056	/* set Outputname (sumpParams.sumpOutfile) *******************************************************/
1057	else if (!strcmp(parmName, "Outputname"))
1058	{
1059	if (expecting == Expecting(EQUALSIGN))
1060	{
1061	expecting = Expecting(ALPHA);
1062	readWord = YES;
1063	}
1064	else if (expecting == Expecting(ALPHA))
1065	{
1066	if(strlen(tkn)>99 && (strchr(tkn,' ')-tkn) > 99 )
1067	{
1068	MrBayesPrint ("%s Maximum allowed length of file name is 99 characters. The given name:\n", spacer);
1069	MrBayesPrint ("%s '%s'\n", spacer,tkn);
1070	return (ERROR);
1071	}
1072	sscanf (tkn, "%s", tempStr);
1073	strcpy (sumpParams.sumpOutfile, tempStr);
1074	MrBayesPrint ("%s Setting sump output file name to \"%s\"\n", spacer, sumpParams.sumpOutfile);
1075	expecting = Expecting(PARAMETER) \| Expecting(SEMICOLON);
1076	}
1077	else
1078	return (ERROR);
1079	}
1080	/* set Relburnin (chainParams.relativeBurnin) ********************************************************/
1081	else if (!strcmp(parmName, "Relburnin"))
1082	{
1083	if (expecting == Expecting(EQUALSIGN))
1084	expecting = Expecting(ALPHA);
1085	else if (expecting == Expecting(ALPHA))
1086	{
1087	if (IsArgValid(tkn, tempStr) == NO_ERROR)
1088	{
1089	if (!strcmp(tempStr, "Yes"))
1090	chainParams.relativeBurnin = YES;
1091	else
1092	chainParams.relativeBurnin = NO;
1093	}
1094	else
1095	{
1096	MrBayesPrint ("%s Invalid argument for Relburnin\n", spacer);
1097	return (ERROR);
1098	}
1099	if (chainParams.relativeBurnin == YES)
1100	MrBayesPrint ("%s Using relative burnin (a fraction of samples discarded).\n", spacer);
1101	else
1102	MrBayesPrint ("%s Using absolute burnin (a fixed number of samples discarded).\n", spacer);
1103	expecting = Expecting(PARAMETER) \| Expecting(SEMICOLON);
1104	}
1105	else
1106	{
1107	return (ERROR);
1108	}
1109	}
1110	/* set Burnin (chainParams.chainBurnIn) ***********************************************************/
1111	else if (!strcmp(parmName, "Burnin"))
1112	{
1113	if (expecting == Expecting(EQUALSIGN))
1114	expecting = Expecting(NUMBER);
1115	else if (expecting == Expecting(NUMBER))
1116	{
1117	sscanf (tkn, "%d", &tempI);
1118	chainParams.chainBurnIn = tempI;
1119	MrBayesPrint ("%s Setting burn-in to %d\n", spacer, chainParams.chainBurnIn);
1120	expecting = Expecting(PARAMETER) \| Expecting(SEMICOLON);
1121	}
1122	else
1123	{
1124	return (ERROR);
1125	}
1126	}
1127	/* set Burninfrac (chainParams.burninFraction) ************************************************************/
1128	else if (!strcmp(parmName, "Burninfrac"))
1129	{
1130	if (expecting == Expecting(EQUALSIGN))
1131	expecting = Expecting(NUMBER);
1132	else if (expecting == Expecting(NUMBER))
1133	{
1134	sscanf (tkn, "%lf", &tempD);
1135	if (tempD < 0.01)
1136	{
1137	MrBayesPrint ("%s Burnin fraction too low (< 0.01)\n", spacer);
1138	return (ERROR);
1139	}
1140	if (tempD > 0.50)
1141	{
1142	MrBayesPrint ("%s Burnin fraction too high (> 0.50)\n", spacer);
1143	return (ERROR);
1144	}
1145	chainParams.burninFraction = tempD;
1146	MrBayesPrint ("%s Setting burnin fraction to %.2f\n", spacer, chainParams.burninFraction);
1147	expecting = Expecting(PARAMETER) \| Expecting(SEMICOLON);
1148	}
1149	else
1150	{
1151	return (ERROR);
1152	}
1153	}
1154	/* set Minprob (sumpParams.minProb) ************************************************************/
1155	else if (!strcmp(parmName, "Minprob"))
1156	{
1157	if (expecting == Expecting(EQUALSIGN))
1158	expecting = Expecting(NUMBER);
1159	else if (expecting == Expecting(NUMBER))
1160	{
1161	sscanf (tkn, "%lf", &tempD);
1162	if (tempD > 0.50)
1163	{
1164	MrBayesPrint ("%s Minprob too high (it should be smaller than 0.50)\n", spacer);
1165	return (ERROR);
1166	}
1167	sumpParams.minProb = tempD;
1168	MrBayesPrint ("%s Setting minprob to %1.3f\n", spacer, sumpParams.minProb);
1169	expecting = Expecting(PARAMETER) \| Expecting(SEMICOLON);
1170	}
1171	else
1172	{
1173	return (ERROR);
1174	}
1175	}
1176	/* set Nruns (sumpParams.numRuns) *******************************************************/
1177	else if (!strcmp(parmName, "Nruns"))
1178	{
1179	if (expecting == Expecting(EQUALSIGN))
1180	expecting = Expecting(NUMBER);
1181	else if (expecting == Expecting(NUMBER))
1182	{
1183	sscanf (tkn, "%d", &tempI);
1184	if (tempI < 1)
1185	{
1186	MrBayesPrint ("%s Nruns must be at least 1\n", spacer);
1187	return (ERROR);
1188	}
1189	else
1190	{
1191	sumpParams.numRuns = tempI;
1192	MrBayesPrint ("%s Setting sump nruns to %d\n", spacer, sumpParams.numRuns);
1193	expecting = Expecting(PARAMETER) \| Expecting(SEMICOLON);
1194	}
1195	}
1196	else
1197	return (ERROR);
1198	}
1199	/* set Hpd (sumpParams.HPD) ********************************************************/
1200	else if (!strcmp(parmName, "Hpd"))
1201	{
1202	if (expecting == Expecting(EQUALSIGN))
1203	expecting = Expecting(ALPHA);
1204	else if (expecting == Expecting(ALPHA))
1205	{
1206	if (IsArgValid(tkn, tempStr) == NO_ERROR)
1207	{
1208	if (!strcmp(tempStr, "Yes"))
1209	sumpParams.HPD = YES;
1210	else
1211	sumpParams.HPD = NO;
1212	}
1213	else
1214	{
1215	MrBayesPrint ("%s Invalid argument for Hpd\n", spacer);
1216	return (ERROR);
1217	}
1218	if (sumpParams.HPD == YES)
1219	MrBayesPrint ("%s Reporting 95 %% region of Highest Posterior Density (HPD).\n", spacer);
1220	else
1221	MrBayesPrint ("%s Reporting median interval containing 95 %% of values.\n", spacer);
1222	expecting = Expecting(PARAMETER) \| Expecting(SEMICOLON);
1223	}
1224	else
1225	{
1226	return (ERROR);
1227	}
1228	}
1229	/* set Allruns (sumpParams.allRuns) ********************************************************/
1230	else if (!strcmp(parmName, "Allruns"))
1231	{
1232	if (expecting == Expecting(EQUALSIGN))
1233	expecting = Expecting(ALPHA);
1234	else if (expecting == Expecting(ALPHA))
1235	{
1236	if (IsArgValid(tkn, tempStr) == NO_ERROR)
1237	{
1238	if (!strcmp(tempStr, "Yes"))
1239	sumpParams.allRuns = YES;
1240	else
1241	sumpParams.allRuns = NO;
1242	}
1243	else
1244	{
1245	MrBayesPrint ("%s Invalid argument for allruns (valid arguments are 'yes' and 'no')\n", spacer);
1246	return (ERROR);
1247	}
1248	if (sumpParams.allRuns == YES)
1249	MrBayesPrint ("%s Setting sump to print information for each run\n", spacer);
1250	else
1251	MrBayesPrint ("%s Setting sump to print only summary information for all runs\n", spacer);
1252	expecting = Expecting(PARAMETER) \| Expecting(SEMICOLON);
1253	}
1254	else
1255	return (ERROR);
1256	}
1257	else
1258	return (ERROR);
1259	}
1260
1261	return (NO_ERROR);
1262
1263	}
1264
1265
1266
1267
1268
1269	int DoSumSsParm (char parmName, char tkn)
1270
1271	{
1272
1273	int tempI;
1274	MrBFlt tempD;
1275	char tempStr[100];
1276
1277	if (expecting == Expecting(PARAMETER))
1278	{
1279	expecting = Expecting(EQUALSIGN);
1280	}
1281	else
1282	{
1283	if (!strcmp(parmName, "Xxxxxxxxxx"))
1284	{
1285	expecting = Expecting(PARAMETER);
1286	expecting \|= Expecting(SEMICOLON);
1287	}
1288	/* set Filename (sumpParams.sumpFileName) ***************************************************/
1289	else if (!strcmp(parmName, "Filename"))
1290	{
1291	if (expecting == Expecting(EQUALSIGN))
1292	{
1293	expecting = Expecting(ALPHA);
1294	readWord = YES;
1295	}
1296	else if (expecting == Expecting(ALPHA))
1297	{
1298	if(strlen(tkn)>99 && (strchr(tkn,' ')-tkn) > 99 )
1299	{
1300	MrBayesPrint ("%s Maximum allowed length of file name is 99 characters. The given name:\n", spacer);
1301	MrBayesPrint ("%s '%s'\n", spacer,tkn);
1302	return (ERROR);
1303	}
1304	sscanf (tkn, "%s", tempStr);
1305	strcpy (sumpParams.sumpFileName, tempStr);
1306	strcpy (sumpParams.sumpOutfile, tempStr);
1307	MrBayesPrint ("%s Setting sump filename and output file name to %s\n", spacer, sumpParams.sumpFileName);
1308	expecting = Expecting(PARAMETER) \| Expecting(SEMICOLON);
1309	}
1310	else
1311	return (ERROR);
1312	}
1313	/* set Outputname (sumpParams.sumpOutfile) *******************************************************/
1314	/*else if (!strcmp(parmName, "Outputname"))
1315	{
1316	if (expecting == Expecting(EQUALSIGN))
1317	{
1318	expecting = Expecting(ALPHA);
1319	readWord = YES;
1320	}
1321	else if (expecting == Expecting(ALPHA))
1322	{
1323	if(strlen(tkn)>99 && (strchr(tkn,' ')-tkn) > 99 )
1324	{
1325	MrBayesPrint ("%s Maximum allowed length of file name is 99 characters. The given name:\n", spacer);
1326	MrBayesPrint ("%s '%s'\n", spacer,tkn);
1327	return (ERROR);
1328	}
1329	sscanf (tkn, "%s", tempStr);
1330	strcpy (sumpParams.sumpOutfile, tempStr);
1331	MrBayesPrint ("%s Setting sump output file name to \"%s\"\n", spacer, sumpParams.sumpOutfile);
1332	expecting = Expecting(PARAMETER) \| Expecting(SEMICOLON);
1333	}
1334	else
1335	return (ERROR);
1336	}*/
1337	/* set Relburnin (chainParams.relativeBurnin) ********************************************************/
1338	else if (!strcmp(parmName, "Relburnin"))
1339	{
1340	if (expecting == Expecting(EQUALSIGN))
1341	expecting = Expecting(ALPHA);
1342	else if (expecting == Expecting(ALPHA))
1343	{
1344	if (IsArgValid(tkn, tempStr) == NO_ERROR)
1345	{
1346	if (!strcmp(tempStr, "Yes"))
1347	chainParams.relativeBurnin = YES;
1348	else
1349	chainParams.relativeBurnin = NO;
1350	}
1351	else
1352	{
1353	MrBayesPrint ("%s Invalid argument for Relburnin\n", spacer);
1354	return (ERROR);
1355	}
1356	if (chainParams.relativeBurnin == YES)
1357	MrBayesPrint ("%s Using relative burnin (a fraction of samples discarded).\n", spacer);
1358	else
1359	MrBayesPrint ("%s Using absolute burnin (a fixed number of samples discarded).\n", spacer);
1360	expecting = Expecting(PARAMETER) \| Expecting(SEMICOLON);
1361	}
1362	else
1363	{
1364	return (ERROR);
1365	}
1366	}
1367	/* set Burnin (chainParams.chainBurnIn) ***********************************************************/
1368	else if (!strcmp(parmName, "Burnin"))
1369	{
1370	if (expecting == Expecting(EQUALSIGN))
1371	expecting = Expecting(NUMBER);
1372	else if (expecting == Expecting(NUMBER))
1373	{
1374	sscanf (tkn, "%d", &tempI);
1375	chainParams.chainBurnIn = tempI;
1376	MrBayesPrint ("%s Setting burn-in to %d\n", spacer, chainParams.chainBurnIn);
1377	expecting = Expecting(PARAMETER) \| Expecting(SEMICOLON);
1378	}
1379	else
1380	{
1381	return (ERROR);
1382	}
1383	}
1384	/* set Burninfrac (chainParams.burninFraction) ************************************************************/
1385	else if (!strcmp(parmName, "Burninfrac"))
1386	{
1387	if (expecting == Expecting(EQUALSIGN))
1388	expecting = Expecting(NUMBER);
1389	else if (expecting == Expecting(NUMBER))
1390	{
1391	sscanf (tkn, "%lf", &tempD);
1392	if (tempD < 0.01)
1393	{
1394	MrBayesPrint ("%s Burnin fraction too low (< 0.01)\n", spacer);
1395	return (ERROR);
1396	}
1397	if (tempD > 0.50)
1398	{
1399	MrBayesPrint ("%s Burnin fraction too high (> 0.50)\n", spacer);
1400	return (ERROR);
1401	}
1402	chainParams.burninFraction = tempD;
1403	MrBayesPrint ("%s Setting burnin fraction to %.2f\n", spacer, chainParams.burninFraction);
1404	expecting = Expecting(PARAMETER) \| Expecting(SEMICOLON);
1405	}
1406	else
1407	{
1408	return (ERROR);
1409	}
1410	}
1411	/* set Nruns (sumssParams.numRuns) *******************************************************/
1412	else if (!strcmp(parmName, "Nruns"))
1413	{
1414	if (expecting == Expecting(EQUALSIGN))
1415	expecting = Expecting(NUMBER);
1416	else if (expecting == Expecting(NUMBER))
1417	{
1418	sscanf (tkn, "%d", &tempI);
1419	if (tempI < 1)
1420	{
1421	MrBayesPrint ("%s Nruns must be at least 1\n", spacer);
1422	return (ERROR);
1423	}
1424	else
1425	{
1426	sumssParams.numRuns = tempI;
1427	MrBayesPrint ("%s Setting sumss nruns to %d\n", spacer, sumssParams.numRuns);
1428	expecting = Expecting(PARAMETER) \| Expecting(SEMICOLON);
1429	}
1430	}
1431	else
1432	return (ERROR);
1433	}
1434	/* set Allruns (sumssParams.allRuns) ********************************************************/
1435	else if (!strcmp(parmName, "Allruns"))
1436	{
1437	if (expecting == Expecting(EQUALSIGN))
1438	expecting = Expecting(ALPHA);
1439	else if (expecting == Expecting(ALPHA))
1440	{
1441	if (IsArgValid(tkn, tempStr) == NO_ERROR)
1442	{
1443	if (!strcmp(tempStr, "Yes"))
1444	sumssParams.allRuns = YES;
1445	else
1446	sumssParams.allRuns = NO;
1447	}
1448	else
1449	{
1450	MrBayesPrint ("%s Invalid argument for allruns (valid arguments are 'yes' and 'no')\n", spacer);
1451	return (ERROR);
1452	}
1453	if (sumssParams.allRuns == YES)
1454	MrBayesPrint ("%s Setting sump to print information for each run\n", spacer);
1455	else
1456	MrBayesPrint ("%s Setting sump to print only summary information for all runs\n", spacer);
1457	expecting = Expecting(PARAMETER) \| Expecting(SEMICOLON);
1458	}
1459	else
1460	return (ERROR);
1461	}
1462	/* set Steptoplot (sumssParams.stepToPlot) *******************************************************/
1463	else if (!strcmp(parmName, "Steptoplot"))
1464	{
1465	if (expecting == Expecting(EQUALSIGN))
1466	expecting = Expecting(NUMBER);
1467	else if (expecting == Expecting(NUMBER))
1468	{
1469	sscanf (tkn, "%d", &tempI);
1470	if (tempI < 0)
1471	{
1472	MrBayesPrint ("%s Steptoplot must be at least 0\n", spacer);
1473	return (ERROR);
1474	}
1475	else
1476	{
1477	sumssParams.stepToPlot = tempI;
1478	MrBayesPrint ("%s Setting sumss steptoplot to %d\n", spacer, sumssParams.stepToPlot);
1479	expecting = Expecting(PARAMETER) \| Expecting(SEMICOLON);
1480	}
1481	}
1482	else
1483	return (ERROR);
1484	}
1485	/* set Smoothing (sumssParams.smoothing ) *******************************************************/
1486	else if (!strcmp(parmName, "Smoothing"))
1487	{
1488	if (expecting == Expecting(EQUALSIGN))
1489	expecting = Expecting(NUMBER);
1490	else if (expecting == Expecting(NUMBER))
1491	{
1492	sscanf (tkn, "%d", &tempI);
1493	if (tempI < 0)
1494	{
1495	MrBayesPrint ("%s Smoothing must be at least 0\n", spacer);
1496	return (ERROR);
1497	}
1498	else
1499	{
1500	sumssParams.smoothing = tempI;
1501	MrBayesPrint ("%s Setting sumss smoothing to %d\n", spacer, sumssParams.smoothing );
1502	expecting = Expecting(PARAMETER) \| Expecting(SEMICOLON);
1503	}
1504	}
1505	else
1506	return (ERROR);
1507	}
1508	/* set Allruns (sumssParams.askForMorePlots) ********************************************************/
1509	else if (!strcmp(parmName, "Askmore"))
1510	{
1511	if (expecting == Expecting(EQUALSIGN))
1512	expecting = Expecting(ALPHA);
1513	else if (expecting == Expecting(ALPHA))
1514	{
1515	if (IsArgValid(tkn, tempStr) == NO_ERROR)
1516	{
1517	if (!strcmp(tempStr, "Yes"))
1518	sumssParams.askForMorePlots = YES;
1519	else
1520	sumssParams.askForMorePlots = NO;
1521	}
1522	else
1523	{
1524	MrBayesPrint ("%s Invalid argument for askmore (valid arguments are 'yes' and 'no')\n", spacer);
1525	return (ERROR);
1526	}
1527	if (sumssParams.askForMorePlots == YES)
1528	MrBayesPrint ("%s Setting sumss to be interactiva by asking for more plots of burn-in or individual steps.\n", spacer);
1529	else
1530	MrBayesPrint ("%s Setting sumss not to be interactive. It will not ask to print more plots.\n", spacer);
1531	expecting = Expecting(PARAMETER) \| Expecting(SEMICOLON);
1532	}
1533	else
1534	return (ERROR);
1535	}
1536	/* set Discardfrac (sumssParams.discardFraction) ************************************************************/
1537	else if (!strcmp(parmName, "Discardfrac"))
1538	{
1539	if (expecting == Expecting(EQUALSIGN))
1540	expecting = Expecting(NUMBER);
1541	else if (expecting == Expecting(NUMBER))
1542	{
1543	sscanf (tkn, "%lf", &tempD);
1544	if (tempD < 0.00)
1545	{
1546	MrBayesPrint ("%s Discard fraction too low (< 0.00)\n", spacer);
1547	return (ERROR);
1548	}
1549	if (tempD > 1.00)
1550	{
1551	MrBayesPrint ("%s Discard fraction too high (> 1.00)\n", spacer);
1552	return (ERROR);
1553	}
1554	sumssParams.discardFraction = tempD;
1555	MrBayesPrint ("%s Setting discard fraction to %.2f\n", spacer, sumssParams.discardFraction);
1556	expecting = Expecting(PARAMETER) \| Expecting(SEMICOLON);
1557	}
1558	else
1559	{
1560	return (ERROR);
1561	}
1562	}
1563	else
1564	return (ERROR);
1565	}
1566
1567	return (NO_ERROR);
1568
1569	}
1570
1571
1572
1573
1574
1575	/* ExamineSumpFile: Collect info on the parameter samples in the file */
1576	int ExamineSumpFile (char fileName, SumpFileInfo fileInfo, char **headerNames, int nHeaders)
1577	{
1578	char sumpTokenP, sumpToken[CMD_STRING_LENGTH], s=NULL, headerLine, t;
1579	int i, lineTerm, inSumpComment, lineNum, lastNonDigitLine, numParamLines, allDigitLine,
1580	lastTokenWasDash, nNumbersOnThisLine, tokenType, burnin, nLines, firstNumCols,
1581	numRows, numColumns;
1582	MrBFlt tempD;
1583	FILE *fp = NULL;
1584
1585
1586	/* open binary file */
1587	if ((fp = OpenBinaryFileR(fileName)) == NULL)
1588	{
1589	/* test for some simple errors */
1590	if (strlen(sumpParams.sumpFileName) > 2)
1591	{
1592	s = sumpParams.sumpFileName + (int) strlen(sumpParams.sumpFileName) - 2;
1593	if (strcmp(s, ".p") == 0)
1594	{
1595	MrBayesPrint ("%s It appears that you need to remove '.p' from the 'Filename' parameter\n", spacer);
1596	MrBayesPrint ("%s Also make sure that 'Nruns' is set correctly\n", spacer);
1597	return ERROR;
1598	}
1599	}
1600	MrBayesPrint ("%s Make sure that 'Nruns' is set correctly\n", spacer);
1601	return ERROR;
1602	}
1603
1604	/* find out what type of line termination is used */
1605	lineTerm = LineTermType (fp);
1606	if (lineTerm != LINETERM_MAC && lineTerm != LINETERM_DOS && lineTerm != LINETERM_UNIX)
1607	{
1608	MrBayesPrint ("%s Unknown line termination\n", spacer);
1609	goto errorExit;
1610	}
1611
1612	/* find length of longest line */
1613	fileInfo->longestLineLength = LongestLine (fp);
1614	fileInfo->longestLineLength += 10; /* better safe than sorry; if you fgets with raw longestLineLength, you run into problems */
1615
1616	/* allocate string long enough to hold a line */
1617	s = (char )SafeMalloc((size_t) (2(fileInfo->longestLineLength + 10) * sizeof(char)));
1618	if (!s)
1619	{
1620	MrBayesPrint ("%s Problem allocating string for reading sump file\n", spacer);
1621	goto errorExit;
1622	}
1623	headerLine = s + fileInfo->longestLineLength + 10;
1624
1625	/* close binary file */
1626	SafeFclose (&fp);
1627
1628	/* open text file */
1629	if ((fp = OpenTextFileR(fileName)) == NULL)
1630	goto errorExit;
1631
1632	/* Check file for appropriate blocks. We want to find the last block
1633	in the file and start from there. */
1634	inSumpComment = NO;
1635	lineNum = lastNonDigitLine = numParamLines = 0;
1636	while (fgets (s, fileInfo->longestLineLength + 2, fp) != NULL)
1637	{
1638	sumpTokenP = &s[0];
1639	allDigitLine = YES;
1640	lastTokenWasDash = NO;
1641	nNumbersOnThisLine = 0;
1642	do {
1643	if(GetToken (sumpToken, &tokenType, &sumpTokenP))
1644	goto errorExit;
1645	/printf ("%s (%d)\n", sumpToken, tokenType);/
1646	if (IsSame("[", sumpToken) == SAME)
1647	inSumpComment = YES;
1648	if (IsSame("]", sumpToken) == SAME)
1649	inSumpComment = NO;
1650
1651	if (inSumpComment == NO)
1652	{
1653	if (tokenType == NUMBER)
1654	{
1655	sscanf (sumpToken, "%lf", &tempD);
1656	if (lastTokenWasDash == YES)
1657	tempD *= -1.0;
1658	nNumbersOnThisLine++;
1659	lastTokenWasDash = NO;
1660	}
1661	else if (tokenType == DASH)
1662	{
1663	lastTokenWasDash = YES;
1664	}
1665	else if (tokenType != UNKNOWN_TOKEN_TYPE)
1666	{
1667	allDigitLine = NO;
1668	lastTokenWasDash = NO;
1669	}
1670	}
1671	} while (*sumpToken);
1672	lineNum++;
1673
1674	if (allDigitLine == NO)
1675	{
1676	lastNonDigitLine = lineNum;
1677	numParamLines = 0;
1678	}
1679	else
1680	{
1681	if (nNumbersOnThisLine > 0)
1682	numParamLines++;
1683	}
1684	}
1685
1686	/* Now, check some aspects of the .p file. */
1687	if (inSumpComment == YES)
1688	{
1689	MrBayesPrint ("%s Unterminated comment in file \"%s\"\n", spacer, fileName);
1690	goto errorExit;
1691	}
1692	if (numParamLines <= 0)
1693	{
1694	MrBayesPrint ("%s No parameters were found in file or there characters not representing a number in last string of the file.\"%s\"\n", spacer, fileName);
1695	goto errorExit;
1696	}
1697
1698	/* calculate burnin */
1699	if ( chainParams.isSS == YES )
1700	{
1701	burnin = 0;
1702	}
1703	else
1704	{
1705	if (chainParams.relativeBurnin == YES)
1706	burnin = (int) (chainParams.burninFraction * numParamLines);
1707	else
1708	burnin = chainParams.chainBurnIn;
1709	}
1710
1711	/* check against burnin */
1712	if (burnin > numParamLines)
1713	{
1714	MrBayesPrint ("%s No parameters can be sampled from file %s as the burnin (%d) exceeds the number of lines in last block (%d)\n",
1715	spacer, fileName, burnin, numParamLines);
1716	MrBayesPrint ("%s Try setting burnin to a number less than %d\n", spacer, numParamLines);
1717	goto errorExit;
1718	}
1719
1720	/* Set some info in fileInfo */
1721	fileInfo->firstParamLine = lastNonDigitLine + burnin;
1722	fileInfo->headerLine = lastNonDigitLine;
1723
1724	/* Calculate and check the number of columns and rows for the file; get header line at the same time */
1725	(void)fseek(fp, 0L, 0);
1726	for (lineNum=0; lineNum<lastNonDigitLine; lineNum++)
1727	if(fgets (s, fileInfo->longestLineLength + 2, fp)==NULL)
1728	goto errorExit;
1729	strcpy(headerLine, s);
1730	for (; lineNum < lastNonDigitLine+burnin; lineNum++)
1731	if(fgets (s, fileInfo->longestLineLength + 2, fp)==NULL)
1732	goto errorExit;
1733
1734	inSumpComment = NO;
1735	nLines = 0;
1736	numRows = numColumns = firstNumCols = 0;
1737	while (fgets (s, fileInfo->longestLineLength + 2, fp) != NULL)
1738	{
1739	sumpTokenP = &s[0];
1740	allDigitLine = YES;
1741	lastTokenWasDash = NO;
1742	nNumbersOnThisLine = 0;
1743	do {
1744	if(GetToken (sumpToken, &tokenType, &sumpTokenP))
1745	goto errorExit;
1746	if (IsSame("[", sumpToken) == SAME)
1747	inSumpComment = YES;
1748	if (IsSame("]", sumpToken) == SAME)
1749	inSumpComment = NO;
1750	if (inSumpComment == NO)
1751	{
1752	if (tokenType == NUMBER)
1753	{
1754	nNumbersOnThisLine++;
1755	lastTokenWasDash = NO;
1756	}
1757	else if (tokenType == DASH)
1758	{
1759	lastTokenWasDash = YES;
1760	}
1761	else if (tokenType != UNKNOWN_TOKEN_TYPE)
1762	{
1763	allDigitLine = NO;
1764	lastTokenWasDash = NO;
1765	}
1766	}
1767	} while (*sumpToken);
1768	lineNum++;
1769	if (allDigitLine == NO)
1770	{
1771	MrBayesPrint ("%s Found a line with non-digit characters (line %d) in file %s\n", spacer, lineNum, fileName);
1772	goto errorExit;
1773	}
1774	else
1775	{
1776	if (nNumbersOnThisLine > 0)
1777	{
1778	nLines++;
1779	if (nLines == 1)
1780	firstNumCols = nNumbersOnThisLine;
1781	else
1782	{
1783	if (nNumbersOnThisLine != firstNumCols)
1784	{
1785	MrBayesPrint ("%s Number of columns is not even (%d in first line and %d in %d line of file %s)\n", spacer, firstNumCols, nNumbersOnThisLine, lineNum, fileName);
1786	goto errorExit;
1787	}
1788	}
1789	}
1790	}
1791	}
1792	fileInfo->numRows = nLines;
1793	fileInfo->numColumns = firstNumCols;
1794
1795	/* set or check headers */
1796	if ((*headerNames) == NULL)
1797	{
1798	GetHeaders (headerNames, headerLine, nHeaders);
1799	if (*nHeaders != fileInfo->numColumns)
1800	{
1801	MrBayesPrint ("%s Expected %d headers but found %d headers\n", spacer, fileInfo->numColumns, *nHeaders);
1802	for (i=0; i<*nHeaders; i++)
1803	SafeFree ((void *) &((headerNames)[i]));
1804	SafeFree ((void *) &(headerNames));
1805	*nHeaders=0;
1806	goto errorExit;
1807	}
1808	}
1809	else
1810	{
1811	if (*nHeaders != fileInfo->numColumns)
1812	{
1813	MrBayesPrint ("%s Expected %d columns but found %d columns\n", spacer, *nHeaders, fileInfo->numColumns);
1814	goto errorExit;
1815	}
1816	for (i=0, t=strtok(headerLine,"\t\n\r"); t!=NULL; t=strtok(NULL,"\t\n\r"), i++)
1817	{
1818	if (i == *nHeaders)
1819	{
1820	MrBayesPrint ("%s Expected %d headers but found more headers.\n",
1821	spacer, fileInfo->numColumns);
1822	goto errorExit;
1823	}
1824	if (strcmp(t,(*headerNames)[i])!=0)
1825	{
1826	MrBayesPrint ("%s Expected header '%s' for column %d but the header for this column was '%s' in file '%s'\n", spacer, (*headerNames)[i], i+1, t, fileName);
1827	MrBayesPrint ("%s It could be that some paramiter values are not numbers and the whole string containing \n",spacer);
1828	MrBayesPrint ("%s this wrongly formated paramiter is treated as a header.\n",spacer);
1829	goto errorExit;
1830	}
1831	}
1832	if (t != NULL)
1833	{
1834	MrBayesPrint ("%s Expected %d headers but found more headers.\n",spacer, fileInfo->numColumns);
1835	goto errorExit;
1836	}
1837	if (i < *nHeaders)
1838	{
1839	MrBayesPrint ("%s Expected header '%s' for column %d but the header for this column was '%s' in file '%s'\n",
1840	spacer, (*headerNames)[i], i+1, t, fileName);
1841	goto errorExit;
1842	}
1843	}
1844
1845	free (s);
1846	fclose(fp);
1847	return (NO_ERROR);
1848
1849	errorExit:
1850
1851	free(s);
1852	fclose(fp);
1853	return (ERROR);
1854	}
1855
1856
1857
1858
1859
1860	/***************************************************
1861	\|
1862	\| FindHeader: Find token in list
1863	\|
1864	----------------------------------------------------*/
1865	int FindHeader (char token, char headerNames, int nHeaders, int index)
1866	{
1867	int i, match=0, nMatches;
1868
1869	*index = -1;
1870	nMatches = 0;
1871	for (i=0; i<nHeaders; i++)
1872	{
1873	if (!strcmp(token,headerNames[i]))
1874	{
1875	nMatches++;
1876	match = i;
1877	}
1878	}
1879
1880	if (nMatches != 1)
1881	return (ERROR);
1882
1883	*index = match;
1884	return (NO_ERROR);
1885	}
1886
1887
1888
1889
1890
1891	/* FreeParameterSamples: Free parameter samples space */
1892	void FreeParameterSamples (ParameterSample *parameterSamples)
1893	{
1894	if (parameterSamples != NULL)
1895	{
1896	free (parameterSamples[0].values[0]);
1897	free (parameterSamples[0].values);
1898	free (parameterSamples);
1899	}
1900	}
1901
1902
1903
1904
1905
1906	/***************************************************
1907	\|
1908	\| GetHeaders: Get headers from headerLine and put
1909	\| them in list while updating nHeaders to reflect
1910	\| the number of headers
1911	\|
1912	----------------------------------------------------*/
1913	int GetHeaders (char **headerNames, char headerLine, int *nHeaders)
1914	{
1915	char *s;
1916
1917	(*nHeaders) = 0;
1918	for (s=strtok(headerLine," \t\n\r"); s!=NULL; s=strtok(NULL," \t\n\r"))
1919	{
1920	if (AddString (headerNames, *nHeaders, s) == ERROR)
1921	{
1922	MrBayesPrint ("%s Error adding header to list of headers \n", spacer, s);
1923	return ERROR;
1924	}
1925	(*nHeaders)++;
1926	}
1927
1928	# if 0
1929	for (i=0; i<(*nHeaders); i++)
1930	printf ("%4d -> '%s'\n", i, headerNames[i]);
1931	# endif
1932
1933	return (NO_ERROR);
1934	}
1935
1936
1937
1938
1939
1940	/* PrintMargLikes: Print marginal likelihoods to screen and to .lstat file */
1941	int PrintMargLikes (char fileName, char headerNames, int nHeaders, ParameterSample parameterSamples, int nRuns, int nSamples)
1942	{
1943	int i, j, len, longestHeader, *sampleCounts=NULL;
1944	char temp[100];
1945	Stat theStats;
1946	FILE *fp;
1947
1948	/* calculate longest header */
1949	longestHeader = 9; /* length of 'parameter' */
1950	for (i=0; i<nHeaders; i++)
1951	{
1952	strcpy (temp, headerNames[i]);
1953	len = (int) strlen(temp);
1954	for (j=0; modelIndicatorParams[j][0]!='\0'; j++)
1955	if (IsSame (temp,modelIndicatorParams[j]) != DIFFERENT)
1956	break;
1957	if (modelIndicatorParams[j][0]!='\0')
1958	continue;
1959	if (!strcmp (temp, "Gen"))
1960	continue;
1961	if (!strcmp (temp, "lnL") == SAME)
1962	continue;
1963	if (len > longestHeader)
1964	longestHeader = len;
1965	}
1966
1967	/* open output file */
1968	strncpy (temp, fileName, 90);
1969	strcat (temp, ".pstat");
1970	fp = OpenNewMBPrintFile (temp);
1971	if (!fp)
1972	return ERROR;
1973
1974	/* print unique identifier to the output file */
1975	if (strlen(stamp) > 1)
1976	fprintf (fp, "[ID: %s]\n", stamp);
1977
1978	/* allocate and set nSamples */
1979	sampleCounts = (int *) SafeCalloc (nRuns, sizeof(int));
1980	if (!sampleCounts)
1981	{
1982	fclose(fp);
1983	return ERROR;
1984	}
1985	for (i=0; i<nRuns; i++)
1986	sampleCounts[i] = nSamples;
1987
1988	/* print the header rows */
1989	MrBayesPrint("\n");
1990	if (sumpParams.HPD == YES)
1991	MrBayesPrint ("%s %*c 95%% HPD Interval\n", spacer, longestHeader, ' ');
1992	else
1993	MrBayesPrint ("%s %*c 95%% Cred. Interval\n", spacer, longestHeader, ' ');
1994	MrBayesPrint ("%s %*c --------------------\n", spacer, longestHeader, ' ');
1995
1996	MrBayesPrint ("%s Parameter%*c Mean Variance Lower Upper Median", spacer, longestHeader-9, ' ');
1997	if (nRuns > 1)
1998	MrBayesPrint (" PSRF+ ");
1999	MrBayesPrint ("\n");
2000
2001	MrBayesPrint ("%s ", spacer);
2002	for (j=0; j<longestHeader+1; j++)
2003	MrBayesPrint ("-");
2004	MrBayesPrint ("-----------------------------------------------------------");
2005	if (nRuns > 1)
2006	MrBayesPrint ("----------");
2007	MrBayesPrint ("\n");
2008	if (nRuns > 1)
2009	MrBayesPrintf (fp, "Parameter\tMean\tVariance\tLower\tUpper\tMedian\tPSRF\n");
2010	else
2011	MrBayesPrintf (fp, "Parameter\tMean\tVariance\tLower\tUpper\tMedian\n");
2012
2013	/* print table values */
2014	for (i=0; i<nHeaders; i++)
2015	{
2016	strcpy (temp, headerNames[i]);
2017	len = (int) strlen(temp);
2018	for (j=0; modelIndicatorParams[j][0]!='\0'; j++)
2019	if (IsSame (temp,modelIndicatorParams[j]) != DIFFERENT)
2020	break;
2021	if (IsSame (temp, "Gen") == SAME)
2022	continue;
2023	if (IsSame (temp, "lnL") == SAME)
2024	continue;
2025
2026	GetSummary (parameterSamples[i].values, nRuns, sampleCounts, &theStats, sumpParams.HPD);
2027
2028	MrBayesPrint ("%s %-*s ", spacer, longestHeader, temp);
2029	MrBayesPrint ("%10.6lf %10.6lf %10.6lf %10.6lf %10.6lf", theStats.mean, theStats.var, theStats.lower, theStats.upper, theStats.median);
2030	MrBayesPrintf (fp, "%s", temp);
2031	MrBayesPrintf (fp, "\t%s", MbPrintNum(theStats.mean));
2032	MrBayesPrintf (fp, "\t%s", MbPrintNum(theStats.var));
2033	MrBayesPrintf (fp, "\t%s", MbPrintNum(theStats.lower));
2034	MrBayesPrintf (fp, "\t%s", MbPrintNum(theStats.upper));
2035	MrBayesPrintf (fp, "\t%s", MbPrintNum(theStats.median));
2036	if (nRuns > 1)
2037	{
2038	if (theStats.PSRF < 0.0)
2039	{
2040	MrBayesPrint (" NA ");
2041	MrBayesPrintf (fp, "NA");
2042	}
2043	else
2044	{
2045	MrBayesPrint (" %7.3lf", theStats.PSRF);
2046	MrBayesPrintf (fp, "\t%s", MbPrintNum(theStats.PSRF));
2047	}
2048	}
2049	MrBayesPrint ("\n");
2050	MrBayesPrintf (fp, "\n");
2051	}
2052	MrBayesPrint ("%s ", spacer);
2053	for (j=0; j<longestHeader+1; j++)
2054	MrBayesPrint ("-");
2055	MrBayesPrint ("-----------------------------------------------------------");
2056	if (nRuns > 1)
2057	MrBayesPrint ("----------");
2058	MrBayesPrint ("\n");
2059	if (nRuns > 1)
2060	{
2061	MrBayesPrint ("%s + Convergence diagnostic (PSRF = Potential Scale Reduction Factor; Gelman\n", spacer);
2062	MrBayesPrint ("%s and Rubin, 1992) should approach 1.0 as runs converge.\n", spacer);
2063	}
2064
2065	fclose (fp);
2066	free (sampleCounts);
2067
2068	return (NO_ERROR);
2069	}
2070
2071
2072
2073
2074
2075	/* PrintModelStats: Print model stats to screen and to .mstat file */
2076	int PrintModelStats (char fileName, char headerNames, int nHeaders, ParameterSample parameterSamples, int nRuns, int nSamples)
2077	{
2078	int i, j, j1, j2, k, longestName, nElements, *modelCounts=NULL;
2079	MrBFlt f, prob=NULL, sum=NULL, ssq=NULL, min=NULL, max=NULL, stddev=NULL;
2080	char temp[100];
2081	FILE *fp;
2082	ModelProb *elem = NULL;
2083
2084	/* nHeaders - is a convenient synonym for number of column headers */
2085
2086	/* check if we have any model indicator variables and also check for longest header */
2087	k = 0;
2088	longestName = 0;
2089	for (i=0; i<nHeaders; i++)
2090	{
2091	for (j=0; strcmp(modelIndicatorParams[j],"")!=0; j++)
2092	{
2093	if (IsSame (headerNames[i], modelIndicatorParams[j]) != DIFFERENT)
2094	{
2095	k++;
2096	for (j1=0; strcmp(modelElementNames[j][j1],"")!=0; j1++)
2097	{
2098	j2 = (int)(strlen(headerNames[i]) + 2 + strlen(modelElementNames[j][j1]));
2099	if (j2 > longestName)
2100	longestName = j2;
2101	}
2102	break;
2103	}
2104	}
2105	}
2106
2107	/* return if nothing to do */
2108	if (k==0)
2109	return NO_ERROR;
2110
2111	/* open output file */
2112	MrBayesPrint ("%s Model probabilities above %1.3lf\n", spacer, sumpParams.minProb);
2113	MrBayesPrint ("%s Estimates saved to file \"%s.mstat\".\n", spacer, sumpParams.sumpOutfile);
2114	strncpy (temp,fileName,90);
2115	strcat (temp, ".mstat");
2116	fp = OpenNewMBPrintFile(temp);
2117	if (!fp)
2118	return ERROR;
2119	MrBayesPrint ("\n");
2120
2121	/* print unique identifier to the output file */
2122	if (strlen(stamp) > 1)
2123	fprintf (fp, "[ID: %s]\n", stamp);
2124
2125	/* print header */
2126	MrBayesPrintf (fp, "\n\n");
2127	if (nRuns == 1)
2128	{
2129	MrBayesPrint ("%s %*c Posterior\n", spacer, longestName-5, ' ');
2130	MrBayesPrint ("%s Model%*c Probability\n", spacer, longestName-5, ' ');
2131	MrBayesPrint ("%s -----", spacer);
2132	for (i=0; i<longestName-5; i++)
2133	MrBayesPrint ("-");
2134	MrBayesPrint ("------------------\n");
2135	MrBayesPrintf (fp, "Model\tProbability\n");
2136	}
2137	else
2138	{
2139	MrBayesPrint ("%s %*c Posterior Standard Min. Max. \n", spacer, longestName-5, ' ');
2140	MrBayesPrint ("%s Model%*c Probability Deviation Probability Probability\n", spacer, longestName-5, ' ');
2141	MrBayesPrint ("%s -----", spacer);
2142	for (i=0; i<longestName-5; i++)
2143	MrBayesPrint ("-");
2144	MrBayesPrint ("---------------------------------------------------------------\n");
2145	MrBayesPrintf (fp, "Model\tProbability\tStd_dev\tMin_prob\tMax_prob\n");
2146	}
2147
2148	/* calculate and print values */
2149	for (i=0; i<nHeaders; i++)
2150	{
2151	for (j=0; modelIndicatorParams[j][0]!='\0'; j++)
2152	if (IsSame (headerNames[i], modelIndicatorParams[j]) != DIFFERENT)
2153	break;
2154	if (modelIndicatorParams[j][0] == '\0')
2155	continue;
2156
2157	for (nElements=0; modelElementNames[j][nElements][0]!='\0'; nElements++)
2158	;
2159
2160	modelCounts = (int *) SafeCalloc (nElements, sizeof(int));
2161	if (!modelCounts)
2162	{
2163	fclose(fp);
2164	return ERROR;
2165	}
2166	prob = (MrBFlt ) SafeCalloc (6nElements, sizeof(MrBFlt));
2167	if (!prob)
2168	{
2169	free (modelCounts);
2170	fclose (fp);
2171	return ERROR;
2172	}
2173	sum = prob + nElements;
2174	ssq = prob + 2*nElements;
2175	stddev = prob + 3*nElements;
2176	min = prob + 4*nElements;
2177	max = prob + 5*nElements;
2178
2179	for (j1=0; j1<nElements; j1++)
2180	min[j1] = 1.0;
2181
2182	for (j1=0; j1<nRuns; j1++)
2183	{
2184	for (j2=0; j2<nElements; j2++)
2185	modelCounts[j2] = 0;
2186	for (j2=0; j2<nSamples; j2++)
2187	modelCounts[(int)(parameterSamples[i].values[j1][j2] + 0.1)]++;
2188	for (j2=0; j2<nElements; j2++)
2189	{
2190	f = (MrBFlt) modelCounts[j2] / (MrBFlt) nSamples;
2191	sum[j2] += f;
2192	ssq[j2] += f*f;
2193	if (f<min[j2])
2194	min[j2] = f;
2195	if (f > max[j2])
2196	max[j2] = f;
2197	}
2198	}
2199
2200	for (j1=0; j1<nElements; j1++)
2201	{
2202	prob[j1] = sum[j1] / (MrBFlt) nRuns;
2203	f = ssq[j1] - (sum[j1] * sum[j1] / (MrBFlt) nRuns);
2204	f /= (nRuns - 1);
2205	if (f <= 0.0)
2206	stddev[j1] = 0.0;
2207	else
2208	stddev[j1] = sqrt (f);
2209	}
2210
2211	elem = (ModelProb *) SafeCalloc (nElements, sizeof(ModelProb));
2212	for (j1=0; j1<nElements; j1++)
2213	{
2214	elem[j1].index = j1;
2215	elem[j1].prob = prob[j1];
2216	}
2217
2218	/* sort in terms of decreasing probabilities */
2219	qsort((void *) elem, (size_t) nElements, (size_t) sizeof(ModelProb), CompareModelProbs);
2220
2221	for (j1=0; j1<nElements; j1++)
2222	{
2223	if (elem[j1].prob <= sumpParams.minProb)
2224	break;
2225
2226	if (nRuns == 1)
2227	{
2228	sprintf (temp, "%s[%s]", headerNames[i], modelElementNames[j][elem[j1].index]);
2229	MrBayesPrint ("%s %-*s %1.3lf\n", spacer, longestName, temp, prob[elem[j1].index]);
2230	MrBayesPrintf (fp, "%s\t%s\n", temp, MbPrintNum(prob[elem[j1].index]));
2231	}
2232	else /* if (nRuns > 1) */
2233	{
2234	sprintf (temp, "%s[%s]", headerNames[i], modelElementNames[j][elem[j1].index]);
2235	MrBayesPrint ("%s %-*s %1.3lf %1.3lf %1.3lf %1.3lf\n",
2236	spacer, longestName, temp, prob[elem[j1].index], stddev[elem[j1].index], min[elem[j1].index], max[elem[j1].index]);
2237	MrBayesPrintf (fp, "%s", temp);
2238	MrBayesPrintf (fp, "\t%s", MbPrintNum(prob[elem[j1].index]));
2239	MrBayesPrintf (fp, "\t%s", MbPrintNum(stddev[elem[j1].index]));
2240	MrBayesPrintf (fp, "\t%s", MbPrintNum(min[elem[j1].index]));
2241	MrBayesPrintf (fp, "\t%s", MbPrintNum(max[elem[j1].index]));
2242	MrBayesPrintf (fp, "\n");
2243	}
2244	}
2245	free(elem);
2246	elem = NULL;
2247	free(modelCounts);
2248	modelCounts = NULL;
2249	free (prob);
2250	prob = NULL;
2251	}
2252
2253	/* print footer */
2254	if (nRuns == 1)
2255	{
2256	MrBayesPrint ("%s -----", spacer);
2257	for (i=0; i<longestName-5; i++)
2258	MrBayesPrint ("-");
2259	MrBayesPrint ("------------------\n\n");
2260	}
2261	else
2262	{
2263	MrBayesPrint ("%s -----", spacer);
2264	for (i=0; i<longestName-5; i++)
2265	MrBayesPrint ("-");
2266	MrBayesPrint ("---------------------------------------------------------------\n\n");
2267	}
2268
2269	/* close output file */
2270	fclose (fp);
2271
2272	return (NO_ERROR);
2273	}
2274
2275
2276
2277
2278
2279	/* PrintOverlayPlot: Print overlay x-y plot of log likelihood vs. generation for several runs */
2280	int PrintOverlayPlot (MrBFlt xVals, MrBFlt yVals, int nRuns, int startingFrom, int nSamples)
2281	{
2282	int i, j, k, k2, n, screenHeight, screenWidth, numY[60], width;
2283	char plotSymbol[15][60];
2284	MrBFlt x, y, minX, maxX, minY, maxY, meanY[60];
2285
2286	if (nRuns == 2)
2287	MrBayesPrint ("\n%s Overlay plot for both runs:\n", spacer);
2288	else
2289	MrBayesPrint ("\n%s Overlay plot for all %d runs:\n", spacer, sumpParams.numRuns);
2290	if (nRuns > 9)
2291	MrBayesPrint ("%s (1 = Run number 1; 2 = Run number 2 etc.; x = Run number 10 or above; * = Several runs)\n", spacer);
2292	else if (nRuns > 2)
2293	MrBayesPrint ("%s (1 = Run number 1; 2 = Run number 2 etc.; * = Several runs)\n", spacer);
2294	else
2295	MrBayesPrint ("%s (1 = Run number 1; 2 = Run number 2; * = Both runs)\n", spacer);
2296
2297	/* print overlay x-y plot of log likelihood vs. generation for all runs */
2298	screenWidth = 60; /* don't change this without changing numY, meanY, and plotSymbol declared above */
2299	screenHeight = 15;
2300
2301	/* find minX, minY, maxX, and maxY over all runs */
2302	minX = minY = 1000000000.0;
2303	maxX = maxY = -1000000000.0;
2304	for (n=0; n<nRuns; n++)
2305	{
2306	for (i=startingFrom; i<startingFrom+nSamples; i++)
2307	{
2308	x = xVals[n][i];
2309	if (x < minX)
2310	minX = x;
2311	if (x > maxX)
2312	maxX = x;
2313	}
2314	}
2315	for (n=0; n<nRuns; n++)
2316	{
2317	y = 0.0;
2318	j = 0;
2319	k2 = 0;
2320	for (i=startingFrom; i<startingFrom+nSamples; i++)
2321	{
2322	x = xVals[n][i];
2323	k = (int) (((x - minX) / (maxX - minX)) * screenWidth);
2324	if (k <= 0)
2325	k = 0;
2326	if (k >= screenWidth)
2327	k = screenWidth - 1;
2328	if (k == j)
2329	{
2330	y += yVals[n][i];
2331	k2 ++;
2332	}
2333	else
2334	{
2335	y /= k2;
2336	if (y < minY)
2337	minY = y;
2338	if (y > maxY)
2339	maxY = y;
2340	k2 = 1;
2341	y = yVals[n][i];
2342	j++;
2343	}
2344	}
2345	if (k2 > 0)
2346	{
2347	y /= k2;
2348	if (y < minY)
2349	minY = y;
2350	if (y > maxY)
2351	maxY = y;
2352	}
2353	}
2354
2355	/* initialize the plot symbols */
2356	for (i=0; i<screenHeight; i++)
2357	for (j=0; j<screenWidth; j++)
2358	plotSymbol[i][j] = ' ';
2359
2360	/* assemble the plot symbols */
2361	for (n=0; n<nRuns; n++)
2362	{
2363	/* find the plot points for this run */
2364	for (i=0; i<screenWidth; i++)
2365	{
2366	numY[i] = 0;
2367	meanY[i] = 0.0;
2368	}
2369	for (i=startingFrom; i<startingFrom+nSamples; i++)
2370	{
2371	x = xVals[n][i];
2372	y = yVals[n][i];
2373	k = (int)(((x - minX) / (maxX - minX)) * screenWidth);
2374	if (k >= screenWidth)
2375	k = screenWidth - 1;
2376	if (k <= 0)
2377	k = 0;
2378	meanY[k] += y;
2379	numY[k]++;
2380	}
2381
2382	/* transfer these points to the overlay */
2383	for (i=0; i<screenWidth; i++)
2384	{
2385	if (numY[i] > 0)
2386	{
2387	k = (int) ((((meanY[i] / numY[i]) - minY)/ (maxY - minY)) * screenHeight);
2388	if (k < 0)
2389	k = 0;
2390	else if (k >= screenHeight)
2391	k = screenHeight - 1;
2392	if (plotSymbol[k][i] == ' ')
2393	{
2394	if (n <= 8)
2395	plotSymbol[k][i] = '1' + n;
2396	else
2397	plotSymbol[k][i] = 'x';
2398	}
2399	else
2400	plotSymbol[k][i] = '*';
2401	}
2402	}
2403	} /* next run */
2404
2405	/* now print the overlay plot */
2406	MrBayesPrint ("\n%s +", spacer);
2407	for (i=0; i<screenWidth; i++)
2408	MrBayesPrint ("-");
2409	MrBayesPrint ("+ %1.2lf\n", maxY);
2410	for (i=screenHeight-1; i>=0; i--)
2411	{
2412	MrBayesPrint ("%s \|", spacer);
2413	for (j=0; j<screenWidth; j++)
2414	{
2415	MrBayesPrint ("%c", plotSymbol[i][j]);
2416	}
2417	MrBayesPrint ("\|\n");
2418	}
2419	MrBayesPrint ("%s +", spacer);
2420	for (i=0; i<screenWidth; i++)
2421	{
2422	if (i % (screenWidth/10) == 0 && i != 0)
2423	MrBayesPrint ("+");
2424	else
2425	MrBayesPrint ("-");
2426	}
2427	MrBayesPrint ("+ %1.2lf\n", minY);
2428	MrBayesPrint ("%s ^", spacer);
2429	for (i=0; i<screenWidth; i++)
2430	MrBayesPrint (" ");
2431	MrBayesPrint ("^\n");
2432	MrBayesPrint ("%s %1.0lf", spacer, minX);
2433	if((int)minX>0)
2434	width=(int)(log10(minX));
2435	else if((int)minX==0)
2436	width=1;
2437	else
2438	width=(int)(log10(-minX))+1;
2439	for (i=0; i<screenWidth-width; i++)
2440	MrBayesPrint (" ");
2441	MrBayesPrint ("%1.0lf\n\n", maxX);
2442
2443	return (NO_ERROR);
2444	}
2445
2446
2447
2448
2449
2450	/* PrintParamStats: Print parameter table (not model indicator params) to screen and .pstat file */
2451	int PrintParamStats (char fileName, char headerNames, int nHeaders, ParameterSample parameterSamples, int nRuns, int nSamples)
2452	{
2453	int i, j, len, longestHeader, *sampleCounts=NULL;
2454	static char *temp=NULL;
2455	char tempf[100];
2456	Stat theStats;
2457	FILE *fp;
2458
2459	/* calculate longest header */
2460	longestHeader = 9; /* length of 'parameter' */
2461	for (i=0; i<nHeaders; i++)
2462	{
2463	SafeStrcpy (&temp, headerNames[i]);
2464	len = (int) strlen(temp);
2465	for (j=0; modelIndicatorParams[j][0]!='\0'; j++)
2466	if (IsSame (temp,modelIndicatorParams[j]) != DIFFERENT)
2467	break;
2468	if (modelIndicatorParams[j][0]!='\0')
2469	continue;
2470	if (!strcmp (temp, "Gen"))
2471	continue;
2472	if (!strcmp (temp, "lnL") == SAME)
2473	continue;
2474	if (len > longestHeader)
2475	longestHeader = len;
2476	}
2477
2478	/* open output file */
2479	strncpy (tempf, fileName, 90);
2480	strcat (tempf, ".pstat");
2481	fp = OpenNewMBPrintFile (tempf);
2482	if (!fp)
2483	return ERROR;
2484
2485	/* print unique identifier to the output file */
2486	if (strlen(stamp) > 1)
2487	fprintf (fp, "[ID: %s]\n", stamp);
2488
2489	/* allocate and set nSamples */
2490	sampleCounts = (int *) SafeCalloc (nRuns, sizeof(int));
2491	if (!sampleCounts)
2492	{
2493	fclose(fp);
2494	return ERROR;
2495	}
2496	for (i=0; i<nRuns; i++)
2497	sampleCounts[i] = nSamples;
2498
2499	/* print the header rows */
2500	MrBayesPrint("\n");
2501	if (sumpParams.HPD == YES)
2502	MrBayesPrint ("%s %*c 95%% HPD Interval\n", spacer, longestHeader, ' ');
2503	else
2504	MrBayesPrint ("%s %*c 95%% Cred. Interval\n", spacer, longestHeader, ' ');
2505	MrBayesPrint ("%s %*c --------------------\n", spacer, longestHeader, ' ');
2506
2507	if (nRuns > 1)
2508	MrBayesPrint ("%s Parameter%c Mean Variance Lower Upper Median min ESS avg ESS PSRF+ ", spacer, longestHeader-9, ' ');
2509	else
2510	MrBayesPrint ("%s Parameter%c Mean Variance Lower Upper Median ESS", spacer, longestHeader-9, ' ');
2511	MrBayesPrint ("\n");
2512
2513	MrBayesPrint ("%s ", spacer);
2514	for (j=0; j<longestHeader+1; j++)
2515	MrBayesPrint ("-");
2516	MrBayesPrint ("---------------------------------------------------------------------");
2517	if (nRuns > 1)
2518	MrBayesPrint ("-------------------");
2519	MrBayesPrint ("\n");
2520	if (nRuns > 1)
2521	MrBayesPrintf (fp, "Parameter\tMean\tVariance\tLower\tUpper\tMedian\tminESS\tavgESS\tPSRF\n");
2522	else
2523	MrBayesPrintf (fp, "Parameter\tMean\tVariance\tLower\tUpper\tMedian\tESS\n");
2524
2525	/* print table values */
2526	for (i=0; i<nHeaders; i++)
2527	{
2528	SafeStrcpy(&temp, headerNames[i]);
2529	len = (int) strlen(temp);
2530	for (j=0; modelIndicatorParams[j][0]!='\0'; j++)
2531	if (IsSame (temp,modelIndicatorParams[j]) != DIFFERENT)
2532	break;
2533	if (modelIndicatorParams[j][0]!='\0')
2534	continue;
2535	if (IsSame (temp, "Gen") == SAME)
2536	continue;
2537	if (IsSame (temp, "lnL") == SAME)
2538	continue;
2539
2540	GetSummary (parameterSamples[i].values, nRuns, sampleCounts, &theStats, sumpParams.HPD);
2541
2542	MrBayesPrint ("%s %-*s ", spacer, longestHeader, temp);
2543	MrBayesPrint ("%10.6lf %10.6lf %10.6lf %10.6lf %10.6lf", theStats.mean, theStats.var, theStats.lower, theStats.upper, theStats.median);
2544	MrBayesPrintf (fp, "%s", temp);
2545	MrBayesPrintf (fp, "\t%s", MbPrintNum(theStats.mean));
2546	MrBayesPrintf (fp, "\t%s", MbPrintNum(theStats.var));
2547	MrBayesPrintf (fp, "\t%s", MbPrintNum(theStats.lower));
2548	MrBayesPrintf (fp, "\t%s", MbPrintNum(theStats.upper));
2549	MrBayesPrintf (fp, "\t%s", MbPrintNum(theStats.median));
2550
2551	if(theStats.minESS == theStats.minESS)
2552	{
2553	MrBayesPrintf (fp, "\t%s", MbPrintNum(theStats.minESS));
2554	MrBayesPrint (" %8.2lf", theStats.minESS);
2555	}
2556	else
2557	{
2558	MrBayesPrint (" NA ");
2559	MrBayesPrintf (fp, "NA");
2560	}
2561	if (nRuns > 1)
2562	{
2563	if(theStats.minESS == theStats.minESS)
2564	{
2565	MrBayesPrint (" %8.2lf", theStats.avrESS);
2566	MrBayesPrintf (fp, "\t%s", MbPrintNum(theStats.avrESS));
2567	}
2568	else
2569	{
2570	MrBayesPrint (" NA ");
2571	MrBayesPrintf (fp, "NA");
2572	}
2573	if (theStats.PSRF < 0.0)
2574	{
2575	MrBayesPrint (" NA ");
2576	MrBayesPrintf (fp, "NA");
2577	}
2578	else
2579	{
2580	MrBayesPrint (" %7.3lf", theStats.PSRF);
2581	MrBayesPrintf (fp, "\t%s", MbPrintNum(theStats.PSRF));
2582	}
2583	}
2584	MrBayesPrint ("\n");
2585	MrBayesPrintf (fp, "\n");
2586	}
2587	MrBayesPrint ("%s ", spacer);
2588	for (j=0; j<longestHeader+1; j++)
2589	MrBayesPrint ("-");
2590	MrBayesPrint ("---------------------------------------------------------------------");
2591	if (nRuns > 1)
2592	MrBayesPrint ("-------------------");
2593	MrBayesPrint ("\n");
2594	if (nRuns > 1)
2595	{
2596	MrBayesPrint ("%s * Convergence diagnostic (ESS = Estimated Sample Size); min and avg values\n", spacer);
2597	MrBayesPrint ("%s correspond to minimal and average ESS among runs. \n", spacer);
2598	MrBayesPrint ("%s ESS value below 100 may indicate that the parameter is undersampled. \n", spacer);
2599	MrBayesPrint ("%s + Convergence diagnostic (PSRF = Potential Scale Reduction Factor; Gelman\n", spacer);
2600	MrBayesPrint ("%s and Rubin, 1992) should approach 1.0 as runs converge.\n", spacer);
2601	}
2602	else
2603	{
2604	MrBayesPrint ("%s * Convergence diagnostic (ESS = Estimated Sample Size); ESS value \n", spacer);
2605	MrBayesPrint ("%s below 100 may indicate that the parameter is undersampled. \n", spacer);
2606	}
2607	MrBayesPrint ("\n\n");
2608
2609	fclose (fp);
2610	free (sampleCounts);
2611	SafeFree ((void **)&temp);
2612
2613	return (NO_ERROR);
2614	}
2615
2616
2617
2618
2619
2620	/* PrintPlot: Print x-y plot of log likelihood vs. generation */
2621	int PrintPlot (MrBFlt xVals, MrBFlt yVals, int numVals)
2622	{
2623	int i, j, k, numY[60], screenWidth, screenHeight;
2624	MrBFlt x, y, minX, maxX, minY, maxY, meanY[60], diff;
2625
2626	/* print x-y plot of log likelihood vs. generation */
2627	screenWidth = 60; /* don't change this without changing numY and meanY, declared above */
2628	screenHeight = 15;
2629
2630	/* find minX and maxX */
2631	minX = xVals[0];
2632	maxX = xVals[0];
2633	for (i=0; i<numVals; i++)
2634	{
2635	x = xVals[i];
2636	if (x < minX)
2637	minX = x;
2638	if (x > maxX)
2639	maxX = x;
2640	}
2641
2642	/* collect Y data */
2643	for (i=0; i<screenWidth; i++)
2644	{
2645	numY[i] = 0;
2646	meanY[i] = 0.0;
2647	}
2648	for (i=0; i<numVals; i++)
2649	{
2650	x = xVals[i];
2651	y = yVals[i];
2652	k = (int)(((x - minX) / (maxX - minX)) * screenWidth);
2653	if (k >= screenWidth)
2654	k = screenWidth - 1;
2655	if (k < 0)
2656	k = 0;
2657	meanY[k] += y;
2658	numY[k]++;
2659	}
2660
2661	/* find minY and maxY */
2662	minY = maxY = meanY[0] / numY[0];
2663	for (i=0; i<screenWidth; i++)
2664	{
2665	if( meanY[i] == 0) /* with some compilers if( NaN < 1 ) is equal true !!! so we realy need this check*/
2666	continue;
2667	meanY[i] /= numY[i];
2668	if (meanY[i] < minY)
2669	minY = meanY[i];
2670	if (meanY[i] > maxY)
2671	maxY = meanY[i];
2672	}
2673
2674	/* find difference */
2675	diff = maxY - minY;
2676
2677	/* print plot */
2678	MrBayesPrint ("\n +");
2679	for (i=0; i<screenWidth; i++)
2680	MrBayesPrint ("-");
2681	MrBayesPrint ("+ %1.3lf\n", maxY);
2682	for (j=screenHeight-1; j>=0; j--)
2683	{
2684	MrBayesPrint (" \|");
2685	for (i=0; i<screenWidth; i++)
2686	{
2687	if (numY[i] > 0)
2688	{
2689	if ((meanY[i] > ((diff/screenHeight)j)+minY && meanY[i] <= ((diff/screenHeight)(j+1))+minY) \|\|
2690	(j == 0 && meanY[i] <= minY))
2691	MrBayesPrint ("*");
2692	else
2693	MrBayesPrint (" ");
2694	}
2695	else
2696	{
2697	MrBayesPrint (" ");
2698	}
2699	}
2700	MrBayesPrint ("\|\n");
2701	}
2702	MrBayesPrint (" +");
2703	for (i=0; i<screenWidth; i++)
2704	{
2705	if (i % (screenWidth/10) == 0 && i != 0)
2706	MrBayesPrint ("+");
2707	else
2708	MrBayesPrint ("-");
2709	}
2710	MrBayesPrint ("+ %1.3lf\n", minY);
2711	MrBayesPrint (" ^");
2712	for (i=0; i<screenWidth; i++)
2713	MrBayesPrint (" ");
2714	MrBayesPrint ("^\n");
2715	MrBayesPrint (" %1.0lf", minX);
2716	if (minX == 0)
2717	j = 1;
2718	else
2719	j = (int)(log10(minX)) + 1;
2720	for (i=0; i<screenWidth-j; i++)
2721	MrBayesPrint (" ");
2722	MrBayesPrint ("%1.0lf\n\n", maxX);
2723
2724	return (NO_ERROR);
2725	}
2726
2727
2728
2729
2730
2731	void PrintPlotHeader (void)
2732	{
2733	MrBayesPrint ("\n");
2734	if (sumpParams.numRuns > 1)
2735	{
2736	MrBayesPrint ("%s Below are rough plots of the generation (x-axis) versus the log \n", spacer);
2737	MrBayesPrint ("%s probability of observing the data (y-axis). You can use these \n", spacer);
2738	MrBayesPrint ("%s graphs to determine what the burn in for your analysis should be. \n", spacer);
2739	MrBayesPrint ("%s When the log probability starts to plateau you may be at station- \n", spacer);
2740	MrBayesPrint ("%s arity. Sample trees and parameters after the log probability \n", spacer);
2741	MrBayesPrint ("%s plateaus. Of course, this is not a guarantee that you are at sta- \n", spacer);
2742	MrBayesPrint ("%s tionarity. Also examine the convergence diagnostics provided by \n", spacer);
2743	MrBayesPrint ("%s the 'sump' and 'sumt' commands for all the parameters in your \n", spacer);
2744	MrBayesPrint ("%s model. Remember that the burn in is the number of samples to dis- \n", spacer);
2745	MrBayesPrint ("%s card. There are a total of ngen / samplefreq samples taken during \n", spacer);
2746	MrBayesPrint ("%s a MCMC analysis. \n", spacer);
2747	}
2748	else
2749	{
2750	MrBayesPrint ("%s Below is a rough plot of the generation (x-axis) versus the log \n", spacer);
2751	MrBayesPrint ("%s probability of observing the data (y-axis). You can use this \n", spacer);
2752	MrBayesPrint ("%s graph to determine what the burn in for your analysis should be. \n", spacer);
2753	MrBayesPrint ("%s When the log probability starts to plateau you may be at station- \n", spacer);
2754	MrBayesPrint ("%s arity. Sample trees and parameters after the log probability \n", spacer);
2755	MrBayesPrint ("%s plateaus. Of course, this is not a guarantee that you are at sta- \n", spacer);
2756	MrBayesPrint ("%s analysis should be. When the log probability starts to plateau \n", spacer);
2757	MrBayesPrint ("%s tionarity. When possible, run multiple analyses starting from dif-\n", spacer);
2758	MrBayesPrint ("%s ferent random trees; if the inferences you make for independent \n", spacer);
2759	MrBayesPrint ("%s analyses are the same, this is reasonable evidence that the chains\n", spacer);
2760	MrBayesPrint ("%s have converged. You can use MrBayes to run several independent \n", spacer);
2761	MrBayesPrint ("%s analyses simultaneously. During such a run, MrBayes will monitor \n", spacer);
2762	MrBayesPrint ("%s the convergence of topologies. After the run has been completed, \n", spacer);
2763	MrBayesPrint ("%s the 'sumt' and 'sump' functions will provide additional conver- \n", spacer);
2764	MrBayesPrint ("%s gence diagnostics for all the parameters in your model. Remember \n", spacer);
2765	MrBayesPrint ("%s that the burn in is the number of samples to discard. There are \n", spacer);
2766	MrBayesPrint ("%s a total of ngen / samplefreq samples taken during a MCMC analysis.\n", spacer);
2767	}
2768	}
2769
2770
2771
2772
2773
2774	/* ReadParamSamples: Read parameter samples from .p file */
2775	int ReadParamSamples (char fileName, SumpFileInfo fileInfo, ParameterSample *parameterSamples, int runNo)
2776	{
2777	char sumpToken[CMD_STRING_LENGTH], s=NULL, p;
2778	int inSumpComment, lineNum, numLinesRead, numLinesToRead, column, lastTokenWasDash,
2779	tokenType;
2780	MrBFlt tempD;
2781	FILE *fp;
2782
2783	/* open file */
2784	if ((fp = OpenTextFileR (fileName)) == NULL)
2785	return ERROR;
2786
2787	/* allocate space for reading lines */
2788	s = (char *) SafeCalloc (fileInfo->longestLineLength + 10, sizeof(char));
2789
2790	/* fast forward to beginning of last unburned parameter line. */
2791	for (lineNum=0; lineNum<fileInfo->firstParamLine; lineNum++)
2792	if (fgets (s, fileInfo->longestLineLength + 5, fp) == 0)
2793	goto errorExit;
2794
2795	/* parse file, line-by-line. We are only parsing lines that have digits that should be read. */
2796	inSumpComment = NO;
2797	numLinesToRead = fileInfo->numRows;
2798	numLinesRead = 0;
2799	while (fgets (s, fileInfo->longestLineLength + 1, fp) != NULL)
2800	{
2801	lastTokenWasDash = NO;
2802	column = 0;
2803	p = s;
2804	do {
2805	if(GetToken (sumpToken, &tokenType, &p))
2806	goto errorExit;
2807	if (IsSame("[", sumpToken) == SAME)
2808	inSumpComment = YES;
2809	if (IsSame("]", sumpToken) == SAME)
2810	inSumpComment = NO;
2811	if (inSumpComment == NO)
2812	{
2813	if (tokenType == NUMBER)
2814	{
2815	/* read the number */
2816	if (column >= fileInfo->numColumns)
2817	{
2818	MrBayesPrint ("%s Too many values read on line %d of file %s\n", spacer, lineNum, fileName);
2819	goto errorExit;
2820	}
2821	sscanf (sumpToken, "%lf", &tempD);
2822	if (lastTokenWasDash == YES)
2823	tempD *= -1.0;
2824	parameterSamples[column].values[runNo][numLinesRead] = tempD;
2825	column++;
2826	lastTokenWasDash = NO;
2827	}
2828	else if (tokenType == DASH)
2829	{
2830	lastTokenWasDash = YES;
2831	}
2832	else if (tokenType != UNKNOWN_TOKEN_TYPE)
2833	{
2834	/* we have a problem */
2835	MrBayesPrint ("%s Line %d of file %s has non-digit characters\n", spacer, lineNum, fileName);
2836	goto errorExit;
2837	}
2838	}
2839	} while (*sumpToken);
2840
2841	lineNum++;
2842	if (column == fileInfo->numColumns)
2843	numLinesRead++;
2844	else if (column != 0)
2845	{
2846	MrBayesPrint ("%s Too few values on line %d of file %s\n", spacer, lineNum, fileName);
2847	goto errorExit;
2848	}
2849	}
2850
2851
2852	/* Check how many parameter line was read in. */
2853	if (numLinesRead != numLinesToRead)
2854	{
2855	MrBayesPrint ("%s Unable to read all lines that should contain parameter samples\n", spacer);
2856	goto errorExit;
2857	}
2858
2859	fclose (fp);
2860	free (s);
2861
2862	return (NO_ERROR);
2863
2864	errorExit:
2865
2866	fclose (fp);
2867	free (s);
2868
2869	return ERROR;
2870	}
2871
2872
2873
2874

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

Context Navigation

source: branches/lib/GDE/MrBAYES/mrbayes_3.2.1/sump.c

Download in other formats: