Context Navigation

← Previous Revision
Next Revision →
Blame
Revision Log

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

Visit:

Last change on this file was 10416, checked in by aboeckma, 12 years ago
Added mr bayes (no menu yet)
File size: 188.3 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 <assert.h>
44
45	#include "mb.h"
46	#include "globals.h"
47	#include "bayes.h"
48	#include "command.h"
49	#include "mbmath.h"
50	#include "mcmc.h"
51	#include "model.h"
52	#include "sump.h"
53	#include "sumt.h"
54	#include "tree.h"
55	#include "utils.h"
56	#if defined(__MWERKS__)
57	#include "SIOUX.h"
58	#endif
59
60	const char* const svnRevisionSumtC="$Rev: 498 $"; /* Revision keyword which is expended/updated by svn on each commit/update*/
61
62	typedef struct partctr
63	{
64	struct partctr left, right;
65	SafeLong *partition;
66	int totCount;
67	int *count;
68	MrBFlt **length;
69	MrBFlt **height;
70	MrBFlt **age;
71	int **nEvents; / nEvents[0,nESets][0,numRuns][0,count[RunID]] */
72	MrBFlt **bRate; / bRate [0,nBSets][0,numRuns][0,count[RunID]] */
73	MrBFlt **bLen; / bLen [0,nBSets][0,numRuns][0,count[RunID]] */
74	MrBFlt *popSize; / popSize[0,numRuns][0,count[RunID]] */
75	}
76	PartCtr;
77
78	typedef struct treectr
79	{
80	struct treectr left, right;
81	int count;
82	int *order;
83	}
84	TreeCtr;
85
86	typedef struct
87	{
88	int longestLineLength;
89	int numTreeBlocks;
90	int lastTreeBlockBegin;
91	int lastTreeBlockEnd;
92	int numTreesInLastBlock;
93	}
94	SumtFileInfo;
95
96	/*#define MAX_PARTITIONS 10000
97	#define MAX_TREES 1000 */
98	#define ALLOC_LEN 100 /* number of values to allocate each time in partition counter nodes */
99
100	#if defined (PRINT_RATEMULTIPLIERS_CPP)
101	FILE *rateMultfp=NULL;
102	#endif
103
104
105	#undef DEBUG_CONTREE
106
107	/* local prototypes */
108	PartCtr AddSumtPartition (PartCtr r, PolyTree t, PolyNode p, int runId);
109	TreeCtr AddSumtTree (TreeCtr r, int *order);
110	PartCtr *AllocPartCtr (void);
111	TreeCtr *AllocTreeCtr (void);
112	void CalculateTreeToTreeDistance (Tree tree1, Tree tree2, MrBFlt d1, MrBFlt d2, MrBFlt *d3);
113	int ConTree (PartCtr **treeParts, int numTreeParts);
114	MrBFlt CppEvolRate (PolyTree t, PolyNode p, int eSet);
115	int ExamineSumtFile (char fileName, SumtFileInfo sumtFileInfo, char treeName, int brlensDef);
116	void FreePartCtr (PartCtr *r);
117	void FreeSumtParams (void);
118	void FreeTreeCtr (TreeCtr *r);
119	int Label (PolyNode p, int addIndex, char label, int maxLength);
120	int OpenBrlensFile (int treeNo);
121	int OpenComptFiles (void);
122	int OpenSumtFiles (int treeNo);
123	void PartCtrUppass (PartCtr r, PartCtr uppass, int index);
124	int PrintBrlensToFile (PartCtr **treeParts, int numTreeParts, int treeNo);
125	void PrintConTree (FILE fp, PolyTree t);
126	void PrintFigTreeConTree (FILE fp, PolyTree t, PartCtr **treeParts);
127	void PrintFigTreeNodeInfo (FILE fp, PartCtr x, MrBFlt length);
128	void PrintSumtTableLine(int numRuns, int rowCount, Stat theStats, MrBFlt numPSRFSamples, MrBFlt maxPSRF, MrBFlt *sumPSRF);
129	void PrintSumtTaxaInfo (void);
130	void Range (MrBFlt vals, int nVals, MrBFlt min, MrBFlt *max);
131	void ResetTaxonSet (void);
132	int ShowConPhylogram (FILE fp, PolyTree t, int screenWidth);
133	void ShowSomeParts (FILE fp, SafeLong p, int offset, int nTaxaToShow);
134	void SortPartCtr (PartCtr **item, int left, int right);
135	void SortTerminalPartCtr (PartCtr **item, int len);
136	void SortTreeCtr (TreeCtr **item, int left, int right);
137	int StoreSumtTree (PackedTree treeList, int index, PolyTree t);
138	void TreeCtrUppass (TreeCtr r, TreeCtr uppass, int index);
139	int TreeProb (void);
140	void WriteConTree (PolyNode p, FILE fp, int showSupport);
141	void WriteFigTreeConTree (PolyNode p, FILE fp, PartCtr **treeParts);
142
143	extern int DoUserTree (void);
144	extern int DoUserTreeParm (char parmName, char tkn);
145	extern int SafeFclose(FILE **);
146
147	extern int inSumtCommand;
148	extern int inComparetreeCommand;
149
150	/* local (to this file) */
151	static int numUniqueSplitsFound, numUniqueTreesFound, numPackedTrees[2], numAsterices; /* length of local to this file variables */
152	static FILE fpParts=NULL, fpTstat=NULL, fpVstat, fpCon=NULL, fpTrees=NULL, fpDists=NULL; /* file pointers */
153	static PartCtr partCtrRoot = NULL; / binary tree for holding splits info */
154	static TreeCtr treeCtrRoot = NULL; / binary tree for holding unique tree info */
155	static PackedTree packedTreeList[2]; / list of trees in packed format */
156
157
158
159
160	PartCtr AddSumtPartition (PartCtr r, PolyTree t, PolyNode p, int runId)
161
162	{
163	int i, n, comp, nLongsNeeded = sumtParams.SafeLongsNeeded;
164
165	if (r == NULL)
166	{
167	/* new partition */
168	/* create a new node */
169	r = AllocPartCtr ();
170	if (r == NULL)
171	return NULL;
172	numUniqueSplitsFound++;
173	for (i=0; i<nLongsNeeded; i++)
174	r->partition[i] = p->partition[i];
175	for (i=0; i<sumtParams.numRuns; i++)
176	r->count[i] = 0;
177	r->left = r->right = NULL;
178	/* record values */
179	if (sumtParams.brlensDef == YES)
180	r->length[runId][0]= p->length;
181	if (sumtParams.isClock == YES)
182	r->height[runId][0]= p->depth;
183	if (sumtParams.isCalibrated == YES)
184	r->age[runId][0]= p->age;
185	for (i=0; i<sumtParams.nESets; i++)
186	r->nEvents[i][runId][0] = t->nEvents[i][p->index];
187	for (i=0; i<sumtParams.nBSets; i++)
188	{
189	r->bLen [i][runId][0] = t->effectiveBrLen[i][p->index];
190	r->bRate[i][runId][0] = t->effectiveBrLen[i][p->index] / p->length;
191	}
192	if (t->popSizeSet == YES)
193	r->popSize[runId][0] = t->popSize[p->index];
194	r->count[runId] ++;
195	r->totCount++;
196	}
197	else
198	{
199	for (i=0; i<nLongsNeeded; i++)
200	{
201	if (r->partition[i] != p->partition[i])
202	break;
203	}
204
205	if (i == nLongsNeeded)
206	comp = 0;
207	else if (r->partition[i] < p->partition[i])
208	comp = -1;
209	else
210	comp = 1;
211
212	if (comp == 0) /* repeated partition */
213	{
214	n = r->count[runId];
215	/* check if we need to allocate more space */
216	if (n % ALLOC_LEN == 0)
217	{
218	/* allocate more space */
219	if (sumtParams.brlensDef == YES)
220	r->length[runId] = (MrBFlt ) SafeRealloc ((void )r->length[runId],(size_t)((n+ALLOC_LEN)*sizeof(MrBFlt)));
221	if (sumtParams.isClock == YES)
222	r->height[runId] = (MrBFlt ) SafeRealloc ((void )r->height[runId],(size_t)((n+ALLOC_LEN)*sizeof(MrBFlt)));
223	if (sumtParams.isCalibrated == YES)
224	r->age[runId] = (MrBFlt ) SafeRealloc ((void )r->age[runId],(size_t)((n+ALLOC_LEN)*sizeof(MrBFlt)));
225	if (sumtParams.nESets > 0)
226	{
227	for (i=0; i<sumtParams.nESets; i++)
228	r->nEvents[i][runId] = (int ) SafeRealloc ((void )r->nEvents[i][runId], (size_t)(n+ALLOC_LEN)*sizeof(int));
229	}
230	if (sumtParams.nBSets > 0)
231	{
232	for (i=0; i<sumtParams.nBSets; i++)
233	{
234	r->bRate[i][runId] = (MrBFlt ) SafeRealloc ((void )r->bRate[i][runId], (size_t)(n+ALLOC_LEN)*sizeof(MrBFlt));
235	r->bLen [i][runId] = (MrBFlt ) SafeRealloc ((void )r->bLen [i][runId], (size_t)(n+ALLOC_LEN)*sizeof(MrBFlt));
236	}
237	}
238	if (sumtParams.popSizeSet == YES)
239	r->popSize[runId] = (MrBFlt ) SafeRealloc ((void )r->popSize[runId],(size_t)((n+ALLOC_LEN)*sizeof(MrBFlt)));
240	}
241	/* record values */
242	r->count[runId]++;
243	r->totCount++;
244	if (sumtParams.brlensDef == YES)
245	r->length[runId][n]= p->length;
246	if (sumtParams.isClock == YES)
247	r->height[runId][n]= p->depth;
248	if (sumtParams.isCalibrated == YES)
249	r->age[runId][n]= p->age;
250	if (sumtParams.nESets > 0)
251	{
252	for (i=0; i<sumtParams.nESets; i++)
253	r->nEvents[i][runId][n] = t->nEvents[i][p->index];
254	}
255	if (sumtParams.nBSets > 0)
256	{
257	for (i=0; i<sumtParams.nBSets; i++)
258	{
259	r->bLen [i][runId][n] = t->effectiveBrLen[i][p->index];
260	r->bRate[i][runId][n] = t->effectiveBrLen[i][p->index] / p->length;
261	}
262	}
263	if (sumtParams.popSizeSet == YES)
264	r->popSize[runId][n] = t->popSize[p->index];
265	}
266	else if (comp < 0) /* greater than -> into left subtree */
267	{
268	if ((r->left = AddSumtPartition (r->left, t, p, runId)) == NULL)
269	{
270	FreePartCtr (r);
271	return NULL;
272	}
273	}
274	else
275	{
276	/* smaller than -> into right subtree */
277	if ((r->right = AddSumtPartition (r->right, t, p, runId)) == NULL)
278	{
279	FreePartCtr (r);
280	return NULL;
281	}
282	}
283	}
284
285	return r;
286	}
287
288
289
290
291
292	TreeCtr AddSumtTree (TreeCtr r, int *order)
293
294	{
295	int i, comp;
296
297	if (r == NULL)
298	{
299	/* new tree */
300	/* create a new node */
301	r = AllocTreeCtr();
302	if (!r)
303	return NULL;
304	numUniqueTreesFound++;
305	for (i=0; i<sumtParams.orderLen; i++)
306	r->order[i] = order[i];
307	r->count = 1;
308	}
309	else
310	{
311	for (i=0; i<sumtParams.orderLen; i++)
312	if (r->order[i] != order[i])
313	break;
314
315	if (i==sumtParams.orderLen)
316	comp = 0;
317	else if (order[i] < r->order[i])
318	comp = 1;
319	else
320	comp = -1;
321
322	if (comp == 0) /* repeated partition */
323	r->count++;
324	else if (comp < 0) /* greater than -> into left subtree */
325	{
326	if ((r->left = AddSumtTree (r->left, order)) == NULL)
327	{
328	FreeTreeCtr (r);
329	return NULL;
330	}
331	}
332	else
333	{
334	/* smaller than -> into right subtree */
335	if ((r->right = AddSumtTree (r->right, order)) == NULL)
336	{
337	FreeTreeCtr (r);
338	return NULL;
339	}
340	}
341	}
342
343	return r;
344	}
345
346
347
348
349
350	/* AllocPartCtr: Allocate space for one partition counter node using info in sumtParams */
351	PartCtr *AllocPartCtr ()
352
353	{
354
355	int i, j;
356	PartCtr *r;
357
358	/* allocate basic stuff */
359	r = (PartCtr *) SafeCalloc ((size_t) 1, sizeof(PartCtr));
360	r->left = r->right = NULL;
361	r->partition = (SafeLong *) SafeCalloc ((size_t) sumtParams.SafeLongsNeeded, sizeof(SafeLong));
362	r->count = (int *) SafeCalloc ((size_t) sumtParams.numRuns, sizeof (int));
363	if (sumtParams.brlensDef)
364	{
365	r->length = (MrBFlt *) SafeCalloc ((size_t) sumtParams.numRuns, sizeof (MrBFlt ));
366	for (i=0; i<sumtParams.numRuns; i++)
367	r->length[i] = (MrBFlt *) SafeCalloc (ALLOC_LEN, sizeof(MrBFlt));
368	}
369	if (sumtParams.isClock)
370	{
371	r->height = (MrBFlt *) SafeCalloc ((size_t) sumtParams.numRuns, sizeof (MrBFlt ));
372	for (i=0; i<sumtParams.numRuns; i++)
373	r->height[i] = (MrBFlt *) SafeCalloc (ALLOC_LEN, sizeof(MrBFlt));
374	r->age = (MrBFlt *) SafeCalloc ((size_t) sumtParams.numRuns, sizeof (MrBFlt ));
375	for (i=0; i<sumtParams.numRuns; i++)
376	r->age[i] = (MrBFlt *) SafeCalloc (ALLOC_LEN, sizeof(MrBFlt));
377	}
378
379	/* allocate relaxed clock parameters: eRate, nEvents, bRate */
380	if (sumtParams.nESets > 0)
381	r->nEvents = (int *) SafeCalloc ((size_t) sumtParams.nESets, sizeof(int ));
382	for (i=0; i<sumtParams.nESets; i++)
383	{
384	r->nEvents[i] = (int *) SafeCalloc ((size_t) sumtParams.numRuns, sizeof(int ));
385	for (j=0; j<sumtParams.numRuns; j++)
386	r->nEvents[i][j] = (int *) SafeCalloc ((size_t) ALLOC_LEN, sizeof(int));
387	}
388	if (sumtParams.nBSets > 0)
389	{
390	r->bLen = (MrBFlt *) SafeCalloc ((size_t) sumtParams.nBSets, sizeof(MrBFlt ));
391	r->bRate = (MrBFlt *) SafeCalloc ((size_t) sumtParams.nBSets, sizeof(MrBFlt ));
392	}
393	for (i=0; i<sumtParams.nBSets; i++)
394	{
395	r->bLen[i] = (MrBFlt *) SafeCalloc ((size_t) sumtParams.numRuns, sizeof(MrBFlt ));
396	r->bRate[i] = (MrBFlt *) SafeCalloc ((size_t) sumtParams.numRuns, sizeof(MrBFlt ));
397	for (j=0; j<sumtParams.numRuns; j++)
398	{
399	r->bLen[i][j] = (MrBFlt *) SafeCalloc ((size_t) ALLOC_LEN, sizeof(MrBFlt));
400	r->bRate[i][j] = (MrBFlt *) SafeCalloc ((size_t) ALLOC_LEN, sizeof(MrBFlt));
401	}
402	}
403	if (sumtParams.popSizeSet == YES)
404	{
405	r->popSize = (MrBFlt *) SafeCalloc ((size_t) sumtParams.numRuns, sizeof (MrBFlt ));
406	for (i=0; i<sumtParams.numRuns; i++)
407	r->popSize[i] = (MrBFlt *) SafeCalloc (ALLOC_LEN, sizeof(MrBFlt));
408	}
409
410	return r;
411	}
412
413
414
415
416
417	/* AllocTreeCtr: Allocate space for a tree counter node using info in sumtParams struct*/
418	TreeCtr *AllocTreeCtr ()
419
420	{
421	TreeCtr *r;
422
423	r = (TreeCtr *) SafeCalloc ((size_t) 1, sizeof(TreeCtr));
424
425	r->left = r->right = NULL;
426
427	r->order = (int *) SafeCalloc ((size_t) sumtParams.orderLen, sizeof(int));
428
429	return r;
430	}
431
432
433
434
435
436	void CalculateTreeToTreeDistance (Tree tree1, Tree tree2, MrBFlt d1, MrBFlt d2, MrBFlt *d3)
437	{
438	int i, j, k;
439	MrBFlt treeLen1=0.0, treeLen2=0.0;
440	TreeNode p, q;
441
442	(d1) = (d2) = (*d3) = 0.0;
443
444	/* set distance-based measures to max value */
445	if (sumtParams.brlensDef == YES)
446	{
447	treeLen1 = TreeLen(tree1);
448	treeLen2 = TreeLen(tree2);
449	(*d2) = treeLen1 + treeLen2;
450	(*d3) = 2.0;
451	}
452
453	/* now we can get distances in a single pass */
454	for (i=0; i<tree1->nNodes; i++)
455	{
456	p = tree1->allDownPass[i];
457	for (j=0; j<tree2->nNodes; j++)
458	{
459	q = tree2->allDownPass[j];
460	for (k=0; k<sumtParams.SafeLongsNeeded; k++)
461	if (p->partition[k] != q->partition[k])
462	break;
463	if (k == sumtParams.SafeLongsNeeded)
464	break;
465	}
466	if (j < tree2->nNodes)
467	{
468	/* match */
469	if (sumtParams.brlensDef == YES)
470	{
471	(*d2) -= (p->length + q->length - fabs(p->length - q->length));
472	(*d3) -= (p->length/treeLen1 + q->length/treeLen2 - fabs(p->length/treeLen1 - q->length/treeLen2));
473	}
474	}
475	else /* if (k < sumtParams.SafeLongsNeeded) */
476	{
477	/* no match */
478	(*d1) += 2.0;
479	}
480	}
481
482	# if 0
483	printf ("DISTANCES: %lf %lf %lf (%lf %lf)\n", d1, d2, *d3, tl1, tl2);
484	for (i=0; i<nnds; i++)
485	{
486	printf ("%4d -- %4d (%lf) %4d (%lf)\n", i, list1[i], lengths1[i], list2[i], lengths2[i]);
487	}
488	# endif
489
490	}
491
492
493
494
495
496	/* ConTree: Construct consensus tree FIXME: numTreeParts is not used*/
497	int ConTree (PartCtr **treeParts, int numTreeParts)
498	{
499	int i, j, targetNode, nBits, isCompat, numTerminalsEncountered;
500	SafeLong x, *partition = NULL;
501	MrBFlt freq, freqInterapted=0;
502	PolyTree t, t2=NULL;
503	PolyNode p, q, r, ql, *pl;
504	PartCtr *part;
505	Stat theStats;
506	int isFirstLoop=1, isInterapted=0;
507
508	/* check that we have at least three species */
509	if (sumtParams.numTaxa < 3)
510	{
511	MrBayesPrint ("%s Too few taxa included to show consensus trees\n", spacer);
512	return ERROR;
513	}
514
515	treeConstruction:
516	/* now, make a consensus tree */
517	/* first allocate and initialize consensus tree */
518	t = AllocatePolyTree(sumtParams.numTaxa);
519	if (!t)
520	{
521	MrBayesPrint ("%s Could not allocate consensus tree\n", spacer);
522	return(ERROR);
523	}
524	t->isRooted = sumtParams.isRooted;
525	t->isClock = sumtParams.isClock;
526	t->isRelaxed = sumtParams.isRelaxed;
527
528	/* initialize consensus tree nodes */
529	for (i=0; i<sumtParams.numTaxa; i++)
530	{
531	t->nodes[i].left = NULL;
532	t->nodes[i].sib = NULL;
533	t->nodes[i].index = i;
534	t->nodes[i].partitionIndex = -1; /* partition ID */
535	t->nodes[i].age = 0.0; /* temporally set to minimum value to allow any insertion in front of the terminal before actual
536	values of age and depth are available */
537	strcpy(t->nodes[i].label, sumtParams.taxaNames[i]);
538	t->nodes[i].depth = 0.0;
539	}
540	for (; i<t->memNodes; i++)
541	{
542	t->nodes[i].left = NULL;
543	t->nodes[i].sib = NULL;
544	t->nodes[i].index = i;
545	t->nodes[i].partitionIndex = -1; /* partition ID */
546	strcpy (t->nodes[i].label, "");
547	}
548
549	/* create bush
550	->x counts number of subtended terminals
551	make sure t->root->left is in outgroup */
552	p = t->root = &t->nodes[sumtParams.numTaxa];
553	p->anc = p->sib = NULL;
554	p->x = sumtParams.numTaxa;
555	p->age = MRBFLT_MAX; /* temporally set to maximum value to allow any insertion in front of the root before actual values of age and depth are available */
556	p->depth = MRBFLT_MAX;
557	j = localOutGroup;
558	q = &t->nodes[j];
559	p->left = q;
560	q->anc = p;
561	q->x = 1;
562	for (i=0; i<sumtParams.numTaxa; i++)
563	{
564	if (i != j)
565	{
566	q->sib = &t->nodes[i];
567	q = q->sib;
568	q->anc = p;
569	q->x = 1;
570	}
571	}
572	q->sib = NULL;
573
574	/* Resolve bush according to partitions.
575	Partitions may include incompatible ones.
576	Partitions must be sorted from most frequent to least frequent
577	for quit test to work when a 50% majority rule tree is requested
578	and in general for consensus tree to be correct. */
579	t->nNodes = sumtParams.numTaxa + 1;
580	t->nIntNodes = 1;
581	if (sumtParams.isRooted == YES)
582	targetNode = 2 * sumtParams.numTaxa - 2;
583	else
584	targetNode = 2 * sumtParams.numTaxa - 3;
585
586	numTerminalsEncountered = 0;
587	for (i=0; i<numUniqueSplitsFound; i++)
588	{
589	/* get partition */
590	part = treeParts[i];
591
592	/* calculate frequency and test if time to quit */
593	if (t->nNodes > targetNode && numTerminalsEncountered == sumtParams.numTaxa)
594	break;
595	freq = (MrBFlt)(part->totCount) / (MrBFlt)(sumtParams.numTreesSampled);
596	if (freq < 0.50 && !strcmp(sumtParams.sumtConType, "Halfcompat"))
597	break;
598
599	/* get partition */
600	partition = part->partition;
601
602	/* count bits in this partition */
603	for (j=nBits=0; j<sumtParams.SafeLongsNeeded; j++)
604	{
605	x = partition[j];
606	for (x = partition[j]; x != 0; x &= (x - 1))
607	nBits++;
608	}
609
610	/* find out if this is an informative partition */
611	if (nBits == sumtParams.numTaxa \|\| nBits == 0)
612	{
613	/* this is the root (for setting age of root node when tree is dated) */
614	q = t->root;
615	q->partitionIndex = i;
616	if (sumtParams.isClock == YES)
617	{
618	GetSummary(part->height, sumtParams.numRuns, part->count, &theStats, sumtParams.HPD);
619	q->depth = theStats.median;
620	for( p = q->left; p!=NULL; p = p->sib )
621	{
622	if( q->depth <= p->depth )
623	break;
624	}
625	assert(p==NULL);/* Root always has 100% freq and it should be older than any other node that has 100% freq. */
626	}
627	if (sumtParams.isCalibrated == YES)
628	{
629	GetSummary(part->age, sumtParams.numRuns, part->count, &theStats, sumtParams.HPD);
630	q->age = theStats.median;
631	for( p = q->left; p!=NULL; p = p->sib )
632	{
633	if( q->age <= p->age )
634	break;
635	}
636	assert(p==NULL);/* Root always has 100% freq and it should be older than any other node that has 100% freq. */
637	}
638	}
639	else if (nBits > 1 && !(nBits == sumtParams.numTaxa - 1 && sumtParams.isRooted == NO))
640	{
641	/* this is an informative partition */
642	/* find anc of partition */
643	j = FirstTaxonInPartition (partition, sumtParams.SafeLongsNeeded);
644	for (p = &t->nodes[j]; p!=NULL; p = p->anc)
645	if (p->x > nBits)
646	break;
647
648	/* do not include if incompatible with ancestor or any of descendants
649	do not check terminals or root because it is
650	redundant and partitions have not necessarily been set for those */
651	isCompat = YES;
652	if (p->anc != NULL && IsPartNested(partition, p->partition, sumtParams.SafeLongsNeeded)==NO)
653	isCompat = NO;
654	else
655	{
656	for (q=p->left; q!=NULL; q=q->sib)
657	{
658	if (q->x > 1 && IsPartCompatible(q->partition, partition, sumtParams.SafeLongsNeeded)==NO)
659	break;
660	}
661	if (q!=NULL)
662	isCompat = NO;
663	}
664
665	if (isCompat == NO)
666	continue;
667
668	/* set new node */
669	q = &t->nodes[t->nNodes];
670	q->partitionIndex = i;
671	q->x = nBits;
672	q->partition = partition;
673	GetSummary(part->length, sumtParams.numRuns, part->count, &theStats, sumtParams.HPD);
674	q->support = freq;
675	q->length = theStats.median;
676	r=NULL;
677	if (sumtParams.isClock == YES)
678	{
679	GetSummary(part->height, sumtParams.numRuns, part->count, &theStats, sumtParams.HPD);
680	q->depth = theStats.median;
681	if( freq < 1.00 )
682	{
683	for( r = p->left; r!=NULL; r = r->sib )
684	{
685	if( IsPartNested(r->partition, partition, sumtParams.SafeLongsNeeded) && r->depth >= q->depth )
686	break; /* child is older then the node we try to add. Not good.*/
687	}
688	if( p->depth <= q->depth )
689	{ /* New node older than the parent. Not good.*/
690	r = p; /* Just to make r!=NULL*/
691	}
692	}
693	}
694	if (sumtParams.isCalibrated == YES)
695	{
696	GetSummary(part->age, sumtParams.numRuns, part->count, &theStats, sumtParams.HPD);
697	q->age = theStats.median;
698	if( freq < 1.00 )
699	{
700	for( r = p->left; r!=NULL; r = r->sib )
701	{
702	if( freq < 1.00 && IsPartNested(r->partition, partition, sumtParams.SafeLongsNeeded) && r->age >= q->age )
703	break; /* child is older then the node we try to add. Not good.*/
704	}
705	if( p->age <= q->age )
706	{ /* New node older than the parent. Not good.*/
707	r = p; /* Just to make r!=NULL*/
708	}
709	}
710	}
711
712	if( r!=NULL && isFirstLoop )
713	{
714	/* cancel the addition of the new node*/
715	isInterapted =1;
716	freqInterapted=freq;
717	break; /* Finish creating the polytree */
718	}
719	t->nNodes++;
720	t->nIntNodes++;
721
722	/* go through descendants of anc */
723	ql = pl = NULL;
724	for (r=p->left; r!=NULL; r=r ->sib)
725	{
726	/* test if r is in the new partition or not */
727	if ((r->x > 1 && IsPartNested(r->partition, partition, sumtParams.SafeLongsNeeded)) \|\| (r->x == 1 && (partition[r->index / nBitsInALong] & (1 << (r->index % nBitsInALong))) != 0))
728	{
729	/* r is in the partition */
730	if (ql == NULL)
731	q->left = r;
732	else
733	ql->sib = r;
734	ql = r;
735	r->anc = q;
736	}
737	else
738	{
739	/* r is not in the partition */
740	if (pl == NULL)
741	p->left = r;
742	else
743	pl->sib = r;
744	pl = r;
745	}
746	}
747	/* terminate new sib-node chain */
748	ql->sib = NULL;
749	/* new node is last in old sib-node chain */
750	pl->sib = q;
751	q->sib = NULL;
752	q->anc = p;
753	}
754	else
755	/* singleton partition */
756	{
757	if (nBits == sumtParams.numTaxa - 1)
758	j = localOutGroup;
759	else
760	j = FirstTaxonInPartition(partition, sumtParams.SafeLongsNeeded); /* nbits == 1 */
761	q = &t->nodes[j];
762	q->partitionIndex = i;
763	q->partition = partition;
764	numTerminalsEncountered++;
765	GetSummary(part->length, sumtParams.numRuns, part->count, &theStats, sumtParams.HPD);
766	q->length = theStats.median;
767	if (sumtParams.isClock == YES)
768	{
769	GetSummary(part->height, sumtParams.numRuns, part->count, &theStats, sumtParams.HPD);
770	q->depth = theStats.median;
771	if(q->anc->depth <= q->depth )
772	{
773	assert(0);/* We never should get here because terminals always have 100% freq and they are younger than any other node that has 100% freq. */
774	}
775	}
776	if (sumtParams.isCalibrated == YES)
777	{
778	GetSummary(part->age, sumtParams.numRuns, part->count, &theStats, sumtParams.HPD);
779	q->age = theStats.median;
780	if(q->anc->age <= q->age )
781	{
782	assert(0);/* We never should get here because terminals always have 100% freq and they are younger than any other node that has 100% freq. */
783	}
784	}
785	}
786	}
787
788	if( isFirstLoop )
789	{
790	t2 = t;
791	if( isInterapted )
792	{
793	isFirstLoop = 0;
794	goto treeConstruction;
795	}
796	}
797
798
799	/* get downpass arrays */
800	GetPolyDownPass(t);
801
802	/* order tips */
803	if (sumtParams.orderTaxa == YES)
804	OrderTips (t);
805
806	if( t!=t2 )
807	{
808	/* get downpass arrays */
809	GetPolyDownPass(t2);
810
811	/* order tips */
812	if (sumtParams.orderTaxa == YES)
813	OrderTips (t2);
814	}
815
816	/* draw tree to stdout and fp */
817	MrBayesPrint ("\n%s Clade credibility values:\n\n", spacer);
818	ShowConTree (stdout, t, 80, YES);
819	if (logToFile == YES)
820	ShowConTree (logFileFp, t, 80, YES);
821	if (sumtParams.brlensDef == YES)
822	{
823	MrBayesPrint ("\n");
824	if (sumtParams.isClock == YES)
825	MrBayesPrint ("%s Phylogram (based on median node depths):\n", spacer);
826	else
827	MrBayesPrint ("%s Phylogram (based on average branch lengths):\n", spacer);
828	if( isInterapted )
829	{
830	MrBayesPrint ("%s Warning. Phylogram containing all nodes with credibility values exceeding\n",spacer);
831	MrBayesPrint ("%s the level set by Contype could not be constructed.\n",spacer);
832	MrBayesPrint ("%s Only nodes with credibility values exceeding %.2f%% (percentage of trees\n", spacer, freqInterapted*100);
833	MrBayesPrint ("%s where the node is present) were included in the phylogram.\n", spacer);
834	}
835	MrBayesPrint ("\n");
836	ShowConPhylogram (stdout, t2, 80);
837	if (logToFile == YES)
838	ShowConPhylogram (logFileFp, t2, 80);
839	}
840
841	/* print taxa block */
842	MrBayesPrintf (fpCon, "begin taxa;\n");
843	MrBayesPrintf (fpCon, "\tdimensions ntax=%d;\n", sumtParams.numTaxa);
844	MrBayesPrintf (fpCon, "\ttaxlabels\n", sumtParams.numTaxa);
845	for (i=0; i<sumtParams.numTaxa; i++)
846	{
847	for (j=0; j<t2->nNodes; j++)
848	if (t2->nodes[j].index == i)
849	break;
850	MrBayesPrintf (fpCon, "\t\t%s\n", t2->nodes[j].label);
851	}
852	MrBayesPrintf (fpCon, "\t\t;\nend;\n");
853
854	MrBayesPrintf (fpCon, "begin trees;\n");
855	MrBayesPrintf (fpCon, "\ttranslate\n");
856	for (i=0; i<sumtParams.numTaxa; i++)
857	{
858	for (j=0; j<t2->nNodes; j++)
859	if (t2->nodes[j].index == i)
860	break;
861	if (i == sumtParams.numTaxa-1)
862	MrBayesPrintf (fpCon, "\t\t%d\t%s\n", t2->nodes[i].index+1, t2->nodes[i].label);
863	else
864	MrBayesPrintf (fpCon, "\t\t%d\t%s,\n", t2->nodes[i].index+1, t2->nodes[i].label);
865	}
866	MrBayesPrintf (fpCon, "\t\t;\n");
867	if (sumtParams.consensusFormat == SIMPLE)
868	PrintConTree(fpCon, t2);
869	else if (sumtParams.consensusFormat == FIGTREE)
870	PrintFigTreeConTree(fpCon, t2, treeParts);
871	MrBayesPrintf (fpCon, "end;\n");
872
873	if( t!=t2 )
874	{
875	FreePolyTree (t2);
876	}
877	/* free memory */
878	FreePolyTree (t);
879
880	return (NO_ERROR);
881	}
882
883
884
885
886
887	MrBFlt CppEvolRate (PolyTree t, PolyNode p, int eSet)
888
889	{
890
891	int i, nEvents;
892	MrBFlt ancRate, branchRate, rate, pos;
893	PolyNode *q;
894
895	nEvents = t->nEvents[eSet][p->index];
896	pos = t->position[eSet][p->index];
897	rate = t->rateMult[eSet][p->index];
898
899	/* note that event positions are from top of branch (more recent, descendant tip) */
900	ancRate = 1.0;
901	// if (t->eType[eSet] == CPPm)
902	// {
903	for (q=p; q->anc != NULL; q=q->anc)
904	{
905	for (i=0; i<t->nEvents[eSet][p->index]; i++)
906	ancRate *= t->rateMult[eSet][p->index][i];
907	}
908	if (nEvents > 0)
909	{
910	branchRate = rate[0] * pos[0];
911	for (i=1; i<nEvents; i++)
912	{
913	branchRate += (pos[i] - pos[i-1]);
914	branchRate *= rate[i];
915	}
916	branchRate += 1.0 - pos[nEvents-1];
917	branchRate *= ancRate;
918	}
919	else
920	branchRate = ancRate;
921	// }
922	/*
923	else if (t->eType[eSet] == CPPi)
924	{
925	for (q=p; q->anc != NULL; q=q->anc)
926	{
927	if (t->nEvents[eSet][p->index]>0)
928	{
929	ancRate = t->rateMult[eSet][p->index][0];
930	break;
931	}
932	}
933	if (nEvents > 0)
934	{
935	branchRate = ancRate * (1.0 - pos[nEvents-1]);
936	for (i=nEvents-2; i>=0; i--)
937	{
938	branchRate += (rate[i+1] * (pos[i+1] - pos[i]));
939	}
940	branchRate += (rate[0] * pos[0]);
941	}
942	else
943	branchRate = ancRate;
944	}
945	*/
946
947	return branchRate;
948	}
949
950
951
952
953
954	int DoCompareTree (void)
955
956	{
957
958	int i, j, k, n, longestLineLength, brlensDef[2], numTreesInLastBlock[2],
959	lastTreeBlockBegin[2], lastTreeBlockEnd[2], xaxis, yaxis, starHolder[80],
960	minNumTrees, screenWidth, screenHeigth, numY[60], nSamples;
961	SafeLong temporarySeed, *mask;
962	PartCtr *x;
963	MrBFlt xProb, yProb, xInc, yInc, xUpper, xLower, yUpper, yLower, dT1=NULL, dT2=NULL, *dT3=NULL, d1, d2, d3,
964	meanY[60], xVal, yVal, minX, minY, maxX, maxY, sums[3];
965	char *s=NULL, prCh, treeName[2][100];
966	FILE *fp;
967	time_t curTime;
968	PartCtr **treeParts=NULL;
969	Tree tree1=NULL, tree2=NULL;
970	SumtFileInfo sumtFileInfo;
971
972	# if defined (MPI_ENABLED)
973	if (proc_id == 0)
974	{
975	# endif
976
977
978	/* Make sure we read trees using DoSumtTree() code instead of with the user tree code */
979	inComparetreeCommand = YES;
980
981	/* set file pointer to NULL */
982	fp = NULL;
983
984	strcpy(treeName[0],"tree"); //in case if paramiter is not specified in a .t file
985	strcpy(treeName[1],"tree");
986
987	/* Check that a data set has been read in. We check taxon names against
988	those read in. */
989	if (isTaxsetDef == NO)
990	{
991	MrBayesPrint ("%s A matrix or set of taxon labels must be specified before comparetree can be used\n", spacer);
992	goto errorExit;
993	}
994
995	/* open output files for summary information (two files); check if we want to overwrite previous results */
996	if (OpenComptFiles () == ERROR)
997	goto errorExit;
998
999	MrBayesPrint ("%s Examining files ...\n", spacer);
1000
1001	/* Examine first file */
1002	if (ExamineSumtFile(comptreeParams.comptFileName1, &sumtFileInfo, treeName[0], &(brlensDef[0])) == ERROR)
1003	return ERROR;
1004
1005	/* Capture info */
1006	longestLineLength = sumtFileInfo.longestLineLength;
1007	numTreesInLastBlock[0] = sumtFileInfo.numTreesInLastBlock;
1008	lastTreeBlockBegin[0] = sumtFileInfo.lastTreeBlockBegin;
1009	lastTreeBlockEnd[0] = sumtFileInfo.lastTreeBlockEnd;
1010
1011	/* Examine second file */
1012	if (ExamineSumtFile(comptreeParams.comptFileName2, &sumtFileInfo, treeName[1], &brlensDef[1]) == ERROR)
1013	return ERROR;
1014
1015	/* Capture info */
1016	if (longestLineLength < sumtFileInfo.longestLineLength)
1017	longestLineLength = sumtFileInfo.longestLineLength;
1018	numTreesInLastBlock[1] = sumtFileInfo.numTreesInLastBlock;
1019	lastTreeBlockBegin[1] = sumtFileInfo.lastTreeBlockBegin;
1020	lastTreeBlockEnd[1] = sumtFileInfo.lastTreeBlockEnd;
1021
1022	/* Check whether we should work with brlens */
1023	if (brlensDef[0] == YES && brlensDef[1] == YES)
1024	sumtParams.brlensDef = YES;
1025	else
1026	sumtParams.brlensDef = NO;
1027
1028	/* Allocate space for command string */
1029	longestLineLength += 10;
1030	s = (char )SafeMalloc((size_t) (longestLineLength sizeof(char)));
1031	if (!s)
1032	{
1033	MrBayesPrint ("%s Problem allocating string for reading tree file\n", spacer);
1034	goto errorExit;
1035	}
1036
1037	/* Allocate space for packed trees */
1038	if (chainParams.relativeBurnin == YES)
1039	{
1040	numPackedTrees[0] = numTreesInLastBlock[0] - (int)(chainParams.burninFraction * numTreesInLastBlock[0]);
1041	numPackedTrees[1] = numTreesInLastBlock[1] - (int)(chainParams.burninFraction * numTreesInLastBlock[1]);
1042	}
1043	else
1044	{
1045	numPackedTrees[0] = numTreesInLastBlock[0] - chainParams.chainBurnIn;
1046	numPackedTrees[1] = numTreesInLastBlock[1] - chainParams.chainBurnIn;
1047	}
1048	if (memAllocs[ALLOC_PACKEDTREES] == YES)
1049	{
1050	MrBayesPrint ("%s packedTreeList is already allocated\n", spacer);
1051	goto errorExit;
1052	}
1053	packedTreeList[0] = (PackedTree *) SafeCalloc(numPackedTrees[0]+numPackedTrees[1], sizeof(PackedTree));
1054	packedTreeList[1] = packedTreeList[0] + numPackedTrees[0];
1055	if (!packedTreeList[0])
1056	{
1057	MrBayesPrint ("%s Problem allocating packed tree list\n", spacer);
1058	goto errorExit;
1059	}
1060	memAllocs[ALLOC_PACKEDTREES] = YES;
1061
1062	/* Tell user we are ready to go */
1063	MrBayesPrint ("%s Summarizing trees in files \"%s\" and \"%s\"\n", spacer,
1064	comptreeParams.comptFileName1,
1065	comptreeParams.comptFileName2);
1066
1067	if (chainParams.relativeBurnin == YES)
1068	MrBayesPrint ("%s Using relative burnin ('relburnin=yes'), discarding the first %.0f %% ('burninfrac=%1.2f') of sampled trees\n",
1069	spacer, chainParams.burninFraction*100.0, chainParams.burninFraction);
1070	else
1071	MrBayesPrint ("%s Using absolute burnin ('relburnin=no'), discarding the first %d ('burnin=%d') sampled trees\n",
1072	spacer, chainParams.chainBurnIn, chainParams.chainBurnIn);
1073
1074	MrBayesPrint ("%s Writing statistics to file %s\n", spacer, comptreeParams.comptOutfile);
1075
1076	/* Set up cheap status bar. */
1077	MrBayesPrint ("\n%s Tree reading status:\n\n", spacer);
1078	MrBayesPrint ("%s 0 10 20 30 40 50 60 70 80 90 100\n", spacer);
1079	MrBayesPrint ("%s v-------v-------v-------v-------v-------v-------v-------v-------v-------v-------v\n", spacer);
1080	MrBayesPrint ("%s *", spacer);
1081	numAsterices = 0;
1082
1083	/* Read file 1 for real */
1084	if ((fp = OpenTextFileR(comptreeParams.comptFileName1)) == NULL)
1085	goto errorExit;
1086
1087	/* ...and fast forward to beginning of last tree block (skipping begin trees). */
1088	for (i=0; i<lastTreeBlockBegin[0] + 1; i++)
1089	{
1090	if( fgets (s, longestLineLength, fp) == NULL )
1091	{
1092	printf("Error in function: %s at line: %d in file: %s", __FUNCTION__, __LINE__, __FILE__);
1093	}
1094	}
1095
1096	/* Calculate burnin */
1097	if (chainParams.relativeBurnin == YES)
1098	comptreeParams.burnin = (int)(chainParams.burninFraction * numTreesInLastBlock[0]);
1099	else
1100	comptreeParams.burnin = chainParams.chainBurnIn;
1101
1102	/* Initialize sumtParams struct */
1103	numUniqueSplitsFound = numUniqueTreesFound = 0;
1104	sumtParams.runId = 0;
1105	strcpy(sumtParams.curFileName, comptreeParams.comptFileName1);
1106	sumtParams.tree = AllocatePolyTree (numTaxa);
1107	AllocatePolyTreePartitions (sumtParams.tree);
1108	sumtParams.numTreesEncountered = sumtParams.numTreesSampled = 0;
1109	sumtParams.numFileTrees = (int ) SafeCalloc (22+2*numTaxa, sizeof(int));
1110	sumtParams.numFileTreesSampled = sumtParams.numFileTrees + sumtParams.numRuns;
1111	sumtParams.order = sumtParams.numFileTrees + 2*sumtParams.numRuns;
1112	sumtParams.absentTaxa = sumtParams.numFileTrees + 2*sumtParams.numRuns + numTaxa;
1113	sumtParams.numTreesInLastBlock = numTreesInLastBlock[0];
1114	if (!sumtParams.numFileTrees)
1115	{
1116	MrBayesPrint ("%s Problems allocating sumtParams.numFileTrees in DoSumt()\n", spacer);
1117	goto errorExit;
1118	}
1119	else
1120	memAllocs[ALLOC_SUMTPARAMS] = YES;
1121
1122	/* ... and parse the file */
1123	expecting = Expecting(COMMAND);
1124	inTreesBlock = YES;
1125	ResetTranslateTable();
1126	for (i=0; i<lastTreeBlockEnd[0] - lastTreeBlockBegin[0] - 1; i++)
1127	{
1128	if( fgets (s, longestLineLength, fp) == NULL )
1129	{
1130	printf("Error in function: %s at line: %d in file: %s", __FUNCTION__, __LINE__, __FILE__);
1131	}
1132	/MrBayesPrint ("%s", s);/
1133	if (ParseCommand (s) == ERROR)
1134	goto errorExit;
1135	}
1136	inTreesBlock = NO;
1137	ResetTranslateTable();
1138
1139	/* Check that at least one tree was read in. */
1140	if (sumtParams.numFileTreesSampled[0] <= 0)
1141	{
1142	MrBayesPrint ("%s No trees read in\n", spacer);
1143	goto errorExit;
1144	}
1145
1146	/* ... and close file */
1147	SafeFclose (&fp);
1148
1149	/* Read file 2 for real */
1150	if ((fp = OpenTextFileR(comptreeParams.comptFileName2)) == NULL)
1151	return ERROR;
1152
1153	/* ...and fast forward to beginning of last tree block. */
1154	for (i=0; i<lastTreeBlockBegin[1] + 1; i++)
1155	{
1156	if( fgets (s, longestLineLength, fp) == NULL )
1157	{
1158	printf("Error in function: %s at line: %d in file: %s", __FUNCTION__, __LINE__, __FILE__);
1159	}
1160	}
1161
1162	/* Renitialize sumtParams struct */
1163	sumtParams.runId = 1;
1164	strcpy (sumtParams.curFileName, comptreeParams.comptFileName2);
1165
1166	/* Calculate burnin */
1167	if (chainParams.relativeBurnin == YES)
1168	comptreeParams.burnin = (int)(chainParams.burninFraction * numTreesInLastBlock[1]);
1169	else
1170	comptreeParams.burnin = chainParams.chainBurnIn;
1171
1172	/* ... and parse the file */
1173	expecting = Expecting(COMMAND);
1174	inTreesBlock = YES;
1175	ResetTranslateTable();
1176	for (i=0; i<lastTreeBlockEnd[1] - lastTreeBlockBegin[1] - 1; i++)
1177	{
1178	if( fgets (s, longestLineLength, fp) == NULL )
1179	{
1180	printf("Error in function: %s at line: %d in file: %s", __FUNCTION__, __LINE__, __FILE__);
1181	}
1182	/MrBayesPrint ("%s", s);/
1183	if (ParseCommand (s) == ERROR)
1184	goto errorExit;
1185	}
1186	inTreesBlock = NO;
1187	ResetTranslateTable();
1188
1189	/* Check that at least one tree was read in. */
1190	if (sumtParams.numFileTreesSampled[1] <= 0)
1191	{
1192	MrBayesPrint ("%s No trees read in\n", spacer);
1193	goto errorExit;
1194	}
1195
1196	/* ... and close file */
1197	SafeFclose (&fp);
1198
1199	/* Now finish cheap status bar. */
1200	if (numAsterices < 80)
1201	for (i=0; i<80 - numAsterices; i++)
1202	MrBayesPrint ("*");
1203	MrBayesPrint ("\n\n");
1204
1205	/* tell user how many trees were successfully read */
1206	/* tell user how many trees were successfully read */
1207	MrBayesPrint ("%s Read %d trees from last tree block of file \"%s\" (sampling %d of them)\n", spacer,
1208	sumtParams.numFileTrees[0],
1209	comptreeParams.comptFileName1,
1210	sumtParams.numFileTreesSampled[0]);
1211	MrBayesPrint ("%s Read %d trees from last tree block of file \"%s\" (sampling %d of them)\n", spacer,
1212	sumtParams.numFileTrees[1],
1213	comptreeParams.comptFileName2,
1214	sumtParams.numFileTreesSampled[1]);
1215
1216	/* Extract partition counter pointers */
1217	treeParts = (PartCtr *) SafeCalloc ((size_t)(numUniqueSplitsFound), sizeof(PartCtr ));
1218	i = 0;
1219	PartCtrUppass(partCtrRoot, treeParts, &i);
1220
1221	/* Sort taxon partitions (clades, splits) ... */
1222	SortPartCtr (treeParts, 0, numUniqueSplitsFound-1);
1223
1224	/* print to screen */
1225	MrBayesPrint (" \n");
1226	MrBayesPrint ("%s General explanation: \n", spacer);
1227	MrBayesPrint (" \n");
1228	MrBayesPrint ("%s In an unrooted tree, a taxon bipartition (split) is specified by removing a \n", spacer);
1229	MrBayesPrint ("%s branch, thereby dividing the species into those to the left and those to the \n", spacer);
1230	MrBayesPrint ("%s right of the branch. Here, taxa to one side of the removed branch are denoted \n", spacer);
1231	MrBayesPrint ("%s '.' and those to the other side are denoted '*'. Specifically, the '.' symbol \n", spacer);
1232	MrBayesPrint ("%s is used for the taxa on the same side as the outgroup. \n", spacer);
1233	MrBayesPrint (" \n");
1234	MrBayesPrint ("%s In a rooted or clock tree, the '*' symbol is simply used to denote the taxa \n", spacer);
1235	MrBayesPrint ("%s that are included in a particular group (clade), that is, all the descendants \n", spacer);
1236	MrBayesPrint ("%s of a particular branch in the tree. Taxa that are not included are denoted \n", spacer);
1237	MrBayesPrint ("%s using the '.' symbol. \n", spacer);
1238	MrBayesPrint (" \n");
1239	MrBayesPrint ("%s The output includes the ID of the encountered clades or splits (sorted from \n", spacer);
1240	MrBayesPrint ("%s highest to lowest probability), the bipartition or clade in '.*' format, \n", spacer);
1241	MrBayesPrint ("%s number of times the bipartition or clade was observed in the first tree file \n", spacer);
1242	MrBayesPrint ("%s served in the first tree file, the number of times the bipartition was, \n", spacer);
1243	MrBayesPrint ("%s observed in the second tree file, the proportion of the time the bipartition \n", spacer);
1244	MrBayesPrint ("%s was found in the first tree file, and the proportion of the time the bi- \n", spacer);
1245	MrBayesPrint ("%s partition was found in the second tree file. \n", spacer);
1246	MrBayesPrint (" \n");
1247	MrBayesPrint ("%s List of taxa in bipartitions: \n", spacer);
1248	MrBayesPrint (" \n");
1249	j = 1;
1250	for (k=0; k<numTaxa; k++)
1251	{
1252	if (sumtParams.absentTaxa[k] == NO && taxaInfo[k].isDeleted == NO)
1253	{
1254	MrBayesPrint ("%s %4d -- %s\n", spacer, j++, taxaNames[k]);
1255	}
1256	}
1257	MrBayesPrint (" \n");
1258	MrBayesPrint ("%s List of taxon bipartitions found in tree files: \n\n", spacer);
1259
1260	i = (int)(log10(sumtParams.numTreesSampled)) - 1;
1261	if (i<1)
1262	i = 1;
1263	j = sumtParams.numTaxa - 8;
1264	if (j < 1)
1265	j = 1;
1266	MrBayesPrint ("%s ID -- Partition%c No1%c No2%*c Freq1 Freq2\n",
1267	spacer, j, ' ', i, ' ', i, ' ');
1268
1269	mask = SafeCalloc (sumtParams.SafeLongsNeeded, sizeof(SafeLong));
1270	for (i=0; i<sumtParams.numTaxa; i++)
1271	SetBit (i, mask);
1272	for (i=0; i<numUniqueSplitsFound; i++)
1273	{
1274	x = treeParts[i];
1275	if (IsBitSet(localOutGroup, x->partition) == YES && sumtParams.isRooted == NO)
1276	FlipBits(x->partition, sumtParams.SafeLongsNeeded, mask);
1277	if ((MrBFlt)x->totCount/(MrBFlt)sumtParams.numTreesSampled >= comptreeParams.minPartFreq)
1278	{
1279	MrBayesPrint ("%s %4d -- ", spacer, i+1);
1280	ShowParts (stdout, x->partition, sumtParams.numTaxa);
1281
1282	j = (int)(log10(sumtParams.numTreesSampled)) + 1;
1283	if (j < 3)
1284	j = 3;
1285	MrBayesPrint (" %d %d %1.3lf %1.3lf\n",
1286	j, x->count[0], j, x->count[1],
1287	(MrBFlt)x->count[0]/(MrBFlt)sumtParams.numFileTreesSampled[0],
1288	(MrBFlt)x->count[1]/(MrBFlt)sumtParams.numFileTreesSampled[1]);
1289
1290	MrBayesPrintf (fpParts, "%d\t%d\t%d\t%1.3lf\t%1.3lf\n",
1291	i+1, x->count[0], x->count[1],
1292	(MrBFlt)x->count[0]/(MrBFlt)sumtParams.numFileTreesSampled[0],
1293	(MrBFlt)x->count[1]/(MrBFlt)sumtParams.numFileTreesSampled[1]);
1294	}
1295	}
1296	free (mask);
1297
1298	/* make a nifty graph plotting frequencies of clades found in the two tree files */
1299	MrBayesPrint (" \n");
1300	MrBayesPrint ("%s Bivariate plot of clade probabilities: \n", spacer);
1301	MrBayesPrint (" \n");
1302	MrBayesPrint ("%s This graph plots the probabilities of clades found in file 1 (the x-axis) \n", spacer);
1303	MrBayesPrint ("%s against the probabilities of the same clades found in file 2 (the y-axis). \n", spacer);
1304	MrBayesPrint (" \n");
1305	xInc = (MrBFlt) (1.0 / 80.0);
1306	yInc = (MrBFlt) (1.0 / 40.0);
1307	yUpper = 1.0;
1308	yLower = yUpper - yInc;
1309	for (yaxis=39; yaxis>=0; yaxis--)
1310	{
1311	xLower = 0.0;
1312	xUpper = xLower + xInc;
1313	for (xaxis=0; xaxis<80; xaxis++)
1314	{
1315	starHolder[xaxis] = 0;
1316	for (i=0; i<numUniqueSplitsFound; i++)
1317	{
1318	x = treeParts[i];
1319	xProb = (MrBFlt)x->count[0]/(MrBFlt)sumtParams.numFileTreesSampled[0];
1320	yProb = (MrBFlt)x->count[1]/(MrBFlt)sumtParams.numFileTreesSampled[1];
1321	if ((xProb > xLower \|\| (xProb == 0.0 && xaxis == 0)) && (xProb <= xUpper \|\| (xProb == 1.0 && xaxis == 79))
1322	&& (yProb > yLower \|\| (yProb == 0.0 && yaxis == 0)) && (yProb <= yUpper \|\| (yProb == 1.0 && yaxis == 39)))
1323	starHolder[xaxis] = 1;
1324	}
1325	xLower += xInc;
1326	xUpper = xLower + xInc;
1327	}
1328
1329	MrBayesPrint ("%s ", spacer);
1330	for (xaxis=0; xaxis<80; xaxis++)
1331	{
1332	prCh = ' ';
1333	if ((xaxis == 0 && yaxis == 0) \|\| (xaxis == 79 && yaxis == 39))
1334	prCh = '+';
1335	else if ((xaxis == 0 && yaxis == 39) \|\| (xaxis == 79 && yaxis == 0))
1336	prCh = '+';
1337	else if ((yaxis == 0 \|\| yaxis == 39) && xaxis > 0 && xaxis < 79)
1338	prCh = '-';
1339	else if ((xaxis == 0 \|\| xaxis == 79) && yaxis > 0 && yaxis < 39)
1340	prCh = '\|';
1341	if (starHolder[xaxis] == 1)
1342	prCh = '*';
1343	MrBayesPrint ("%c", prCh);
1344	}
1345	if (yaxis == 39)
1346	MrBayesPrint (" 1.00\n");
1347	else if (yaxis == 0)
1348	MrBayesPrint (" 0.00\n");
1349	else
1350	MrBayesPrint ("\n");
1351
1352	yUpper -= yInc;
1353	yLower = yUpper - yInc;
1354	}
1355
1356	MrBayesPrint ("%s ^ ^\n", spacer);
1357	MrBayesPrint ("%s 0.00 1.00\n", spacer);
1358
1359	/* get tree-to-tree distances: first allocate some space */
1360	minNumTrees = sumtParams.numFileTreesSampled[0];
1361	if (sumtParams.numFileTreesSampled[1] < minNumTrees)
1362	minNumTrees = sumtParams.numFileTreesSampled[1];
1363	dT1 = (MrBFlt )SafeMalloc((size_t) (3 minNumTrees * sizeof(MrBFlt)));
1364	tree1 = AllocateFixedTree (sumtParams.numTaxa, sumtParams.isRooted);
1365	tree2 = AllocateFixedTree (sumtParams.numTaxa, sumtParams.isRooted);
1366	if (!dT1 \|\| !tree1 \|\| !tree2)
1367	{
1368	MrBayesPrint ("%s Problem allocating topological distances\n", spacer);
1369	goto errorExit;
1370	}
1371	dT2 = dT1 + minNumTrees;
1372	dT3 = dT2 + minNumTrees;
1373
1374	for (i=0; i<minNumTrees; i++)
1375	{
1376	if (sumtParams.isRooted == NO)
1377	{
1378	RetrieveUTree (tree1, packedTreeList[0][i].order, packedTreeList[0][i].brlens);
1379	RetrieveUTree (tree2, packedTreeList[1][i].order, packedTreeList[1][i].brlens);
1380	}
1381	else
1382	{
1383	RetrieveRTree (tree1, packedTreeList[0][i].order, packedTreeList[0][i].brlens);
1384	RetrieveRTree (tree2, packedTreeList[1][i].order, packedTreeList[1][i].brlens);
1385	}
1386	/* Allocate and set partitions now that we have a tree */
1387	if (i == 0)
1388	{
1389	AllocateTreePartitions(tree1);
1390	AllocateTreePartitions(tree2);
1391	}
1392	else
1393	{
1394	ResetTreePartitions(tree1);
1395	ResetTreePartitions(tree2);
1396	}
1397	CalculateTreeToTreeDistance (tree1, tree2, &d1, &d2, &d3);
1398	dT1[i] = d1;
1399	dT2[i] = d2;
1400	dT3[i] = d3;
1401	}
1402
1403	for (i=0; i<minNumTrees; i++)
1404	{
1405	/MrBayesPrint ("%s %4d -- %lf %lf %lf\n", spacer, i+1, dT1[i], dT2[i], dT3[i]);/
1406	if (sumtParams.brlensDef == YES)
1407	MrBayesPrintf (fpDists, "%d\t%lf\t%lf\t%lf\n", i+1, dT1[i], dT2[i], dT3[i]);
1408	else
1409	MrBayesPrintf (fpDists, "%d\t%lf\n", i+1, dT1[i]);
1410	}
1411
1412	/* print x-y plot of log likelihood vs. generation */
1413	MrBayesPrint ("\n");
1414	MrBayesPrint ("%s Rough plots of generation (x-axis) versus the measures of tree- \n", spacer);
1415	MrBayesPrint ("%s to-tree distances (y-axis). \n", spacer);
1416	MrBayesPrint ("\n");
1417	MrBayesPrint ("%s Distance(Robinson-Foulds):\n", spacer);
1418	MrBayesPrint ("\n");
1419	screenWidth = 60; /* don't change this without changing numY and meanY, declared above */
1420	screenHeigth = 15;
1421	minX = minY = 1000000000.0;
1422	maxX = maxY = -1000000000.0;
1423	for (i=0; i<minNumTrees; i++)
1424	{
1425	xVal = (MrBFlt) (i + chainParams.chainBurnIn);
1426	yVal = dT1[i];
1427	if (xVal < minX)
1428	minX = xVal;
1429	if (yVal < minY)
1430	minY = yVal;
1431	if (xVal > maxX)
1432	maxX = xVal;
1433	if (yVal > maxY)
1434	maxY = yVal;
1435	}
1436	for (i=0; i<screenWidth; i++)
1437	{
1438	numY[i] = 0;
1439	meanY[i] = 0.0;
1440	}
1441	for (i=0; i<minNumTrees; i++)
1442	{
1443	xVal = (MrBFlt) (i + chainParams.chainBurnIn);
1444	yVal = dT1[i];
1445	k = (int)(((xVal - minX) / (maxX - minX)) * screenWidth);
1446	if (k >= screenWidth)
1447	k = screenWidth - 1;
1448	meanY[k] += yVal;
1449	numY[k]++;
1450	}
1451	MrBayesPrint ("\n%s +", spacer);
1452	for (i=0; i<screenWidth; i++)
1453	MrBayesPrint ("-");
1454	MrBayesPrint ("+ %1.2lf\n", maxY);
1455	for (j=screenHeigth-1; j>=0; j--)
1456	{
1457	MrBayesPrint ("%s \|", spacer);
1458	for (i=0; i<screenWidth; i++)
1459	{
1460	if (numY[i] > 0)
1461	{
1462	if (meanY[i] / numY[i] > (((maxY - minY)/screenHeigth)j)+minY && meanY[i] / numY[i] <= (((maxY - minY)/screenHeigth)(j+1))+minY)
1463	MrBayesPrint ("*");
1464	else
1465	MrBayesPrint (" ");
1466	}
1467	else
1468	{
1469	MrBayesPrint (" ");
1470	}
1471	}
1472	MrBayesPrint ("\|\n");
1473	}
1474	MrBayesPrint ("%s +", spacer);
1475	for (i=0; i<screenWidth; i++)
1476	{
1477	if (numY[i] > 0 && meanY[i] / numY[i] <= minY)
1478	MrBayesPrint ("*");
1479	else if (i % (screenWidth/10) == 0 && i != 0)
1480	MrBayesPrint ("+");
1481	else
1482	MrBayesPrint ("-");
1483	}
1484	MrBayesPrint ("+ %1.2lf\n", minY);
1485	MrBayesPrint ("%s ^", spacer);
1486	for (i=0; i<screenWidth; i++)
1487	MrBayesPrint (" ");
1488	MrBayesPrint ("^\n");
1489	MrBayesPrint ("%s %1.0lf", spacer, minX);
1490	for (i=0; i<screenWidth; i++)
1491	MrBayesPrint (" ");
1492	MrBayesPrint ("%1.0lf\n\n", maxX);
1493
1494	if (sumtParams.brlensDef == YES)
1495	{
1496	for (n=0; n<2; n++)
1497	{
1498	MrBayesPrint ("\n");
1499	if (n == 0)
1500	MrBayesPrint ("%s Distance(Robinson-Foulds with branch lengths):\n", spacer);
1501	else
1502	MrBayesPrint ("%s Distance(Robinson-Foulds with scaled branch lengths):\n", spacer);
1503	MrBayesPrint ("\n");
1504	screenWidth = 60; /* don't change this without changing numY and meanY, declared above */
1505	screenHeigth = 15;
1506	minX = minY = 1000000000.0;
1507	maxX = maxY = -1000000000.0;
1508	for (i=0; i<minNumTrees; i++)
1509	{
1510	xVal = (MrBFlt) (i + comptreeParams.burnin);
1511	if (n == 0)
1512	yVal = dT2[i];
1513	else
1514	yVal = dT3[i];
1515	if (xVal < minX)
1516	minX = xVal;
1517	if (yVal < minY)
1518	minY = yVal;
1519	if (xVal > maxX)
1520	maxX = xVal;
1521	if (yVal > maxY)
1522	maxY = yVal;
1523	}
1524	for (i=0; i<screenWidth; i++)
1525	{
1526	numY[i] = 0;
1527	meanY[i] = 0.0;
1528	}
1529	for (i=0; i<minNumTrees; i++)
1530	{
1531	xVal = (MrBFlt) (i + comptreeParams.burnin);
1532	if (n == 0)
1533	yVal = dT2[i];
1534	else
1535	yVal = dT3[i];
1536	k = (int)(((xVal - minX) / (maxX - minX)) * screenWidth);
1537	if (k >= screenWidth)
1538	k = screenWidth - 1;
1539	meanY[k] += yVal;
1540	numY[k]++;
1541	}
1542	MrBayesPrint ("\n%s +", spacer);
1543	for (i=0; i<screenWidth; i++)
1544	MrBayesPrint ("-");
1545	MrBayesPrint ("+ %1.2lf\n", maxY);
1546	for (j=screenHeigth-1; j>=0; j--)
1547	{
1548	MrBayesPrint ("%s \|", spacer);
1549	for (i=0; i<screenWidth; i++)
1550	{
1551	if (numY[i] > 0)
1552	{
1553	if (meanY[i] / numY[i] > (((maxY - minY)/screenHeigth)j)+minY && meanY[i] / numY[i] <= (((maxY - minY)/screenHeigth)(j+1))+minY)
1554	MrBayesPrint ("*");
1555	else
1556	MrBayesPrint (" ");
1557	}
1558	else
1559	{
1560	MrBayesPrint (" ");
1561	}
1562	}
1563	MrBayesPrint ("\|\n");
1564	}
1565	MrBayesPrint ("%s +", spacer);
1566	for (i=0; i<screenWidth; i++)
1567	{
1568	if (numY[i] > 0 && meanY[i] / numY[i] <= minY)
1569	MrBayesPrint ("*");
1570	else if (i % (screenWidth/10) == 0 && i != 0)
1571	MrBayesPrint ("+");
1572	else
1573	MrBayesPrint ("-");
1574	}
1575	MrBayesPrint ("+ %1.2lf\n", minY);
1576	MrBayesPrint ("%s ^", spacer);
1577	for (i=0; i<screenWidth; i++)
1578	MrBayesPrint (" ");
1579	MrBayesPrint ("^\n");
1580	MrBayesPrint ("%s %1.0lf", spacer, minX);
1581	for (i=0; i<screenWidth; i++)
1582	MrBayesPrint (" ");
1583	MrBayesPrint ("%1.0lf\n\n", maxX);
1584	}
1585	}
1586
1587	/* calculate average tree-to-tree distances */
1588	curTime = time(NULL);
1589	temporarySeed = (SafeLong)curTime;
1590	if (temporarySeed < 0)
1591	temporarySeed = -temporarySeed;
1592	sums[0] = sums[1] = sums[2] = 0.0;
1593	nSamples = 1000;
1594	for (n=0; n<nSamples; n++)
1595	{
1596	i = (int) RandomNumber(&temporarySeed) * minNumTrees;
1597	j = (int) RandomNumber(&temporarySeed) * minNumTrees;
1598	if (sumtParams.isRooted == NO)
1599	{
1600	RetrieveUTree (tree1, packedTreeList[0][i].order, packedTreeList[0][i].brlens);
1601	RetrieveUTree (tree2, packedTreeList[1][j].order, packedTreeList[1][j].brlens);
1602	}
1603	else /* if (sumtParams.isRooted == YES) */
1604	{
1605	RetrieveRTree (tree1, packedTreeList[0][i].order, packedTreeList[0][i].brlens);
1606	RetrieveRTree (tree2, packedTreeList[1][j].order, packedTreeList[1][j].brlens);
1607	}
1608	CalculateTreeToTreeDistance (tree1, tree2, &d1, &d2, &d3);
1609	sums[0] += d1;
1610	sums[1] += d2;
1611	sums[2] += d3;
1612	}
1613	MrBayesPrint ("%s Mean tree-to-tree distances, based on %d trees randomly sampled from both files:\n\n", spacer, nSamples);
1614	MrBayesPrint ("%s Mean(Robinson-Foulds) = %1.3lf\n", spacer, sums[0]/nSamples);
1615	if (sumtParams.brlensDef == YES)
1616	{
1617	MrBayesPrint ("%s Mean(Robinson-Foulds with branch lengths) = %1.3lf\n", spacer, sums[1]/nSamples);
1618	MrBayesPrint ("%s Mean(Robinson-Foulds with scaled branch lengths) = %1.3lf\n", spacer, sums[2]/nSamples);
1619	}
1620	MrBayesPrint ("\n");
1621
1622	/* free memory and file pointers */
1623	free(s);
1624	free (dT1);
1625	FreeTree (tree1);
1626	FreeTree (tree2);
1627	if (memAllocs[ALLOC_PACKEDTREES] == YES)
1628	{
1629	for (i=0; i<numPackedTrees[0]+numPackedTrees[1]; i++)
1630	{
1631	free(packedTreeList[0][i].order);
1632	free(packedTreeList[0][i].brlens);
1633	}
1634	free(packedTreeList[0]);
1635	memAllocs[ALLOC_PACKEDTREES] = NO;
1636	}
1637
1638	/* free sumtParams */
1639	FreeSumtParams();
1640
1641	/* close files */
1642	SafeFclose (&fp);
1643	SafeFclose (&fpParts);
1644	SafeFclose (&fpDists);
1645
1646	/* free pointer array to partitions, part and tree counters */
1647	free (treeParts);
1648	FreePartCtr (partCtrRoot);
1649	FreeTreeCtr (treeCtrRoot);
1650	partCtrRoot = NULL;
1651	treeCtrRoot = NULL;
1652
1653	/* reset taxon set */
1654	ResetTaxonSet();
1655
1656	# if defined (MPI_ENABLED)
1657	}
1658	# endif
1659
1660	expecting = Expecting(COMMAND);
1661	inComparetreeCommand = NO;
1662	return (NO_ERROR);
1663
1664	/* error exit */
1665	errorExit:
1666	if (s) free(s);
1667
1668	/* free sumtParams */
1669	FreeSumtParams();
1670
1671	free (dT1);
1672	FreeTree (tree1);
1673	FreeTree (tree2);
1674	if (memAllocs[ALLOC_PACKEDTREES] == YES)
1675	{
1676	for (i=0; i<numPackedTrees[0]+numPackedTrees[1]; i++)
1677	{
1678	free(packedTreeList[0][i].order);
1679	free(packedTreeList[0][i].brlens);
1680	}
1681	free(packedTreeList[0]);
1682	memAllocs[ALLOC_PACKEDTREES] = NO;
1683	}
1684
1685	/* reset taxon set */
1686	ResetTaxonSet();
1687
1688	/* free pointer array to partitions, part and tree counters */
1689	free (treeParts);
1690	FreePartCtr (partCtrRoot);
1691	FreeTreeCtr (treeCtrRoot);
1692	partCtrRoot = NULL;
1693	treeCtrRoot = NULL;
1694
1695	SafeFclose (&fp);
1696	SafeFclose (&fpParts);
1697	SafeFclose (&fpDists);
1698	strcpy (spacer, "");
1699
1700	expecting = Expecting(COMMAND);
1701	inComparetreeCommand = NO;
1702
1703	return (ERROR);
1704
1705	}
1706
1707
1708
1709
1710
1711	int DoCompareTreeParm (char parmName, char tkn)
1712
1713	{
1714
1715	int tempI;
1716	MrBFlt tempD;
1717	char tempStr[100];
1718
1719	if (expecting == Expecting(PARAMETER))
1720	{
1721	expecting = Expecting(EQUALSIGN);
1722	}
1723	else
1724	{
1725	if (!strcmp(parmName, "Xxxxxxxxxx"))
1726	{
1727	expecting = Expecting(PARAMETER);
1728	expecting \|= Expecting(SEMICOLON);
1729	}
1730	/* set Filename (comptreeParams.comptFileName1) *************************************************/
1731	else if (!strcmp(parmName, "Filename1"))
1732	{
1733	if (expecting == Expecting(EQUALSIGN))
1734	{
1735	expecting = Expecting(ALPHA);
1736	readWord = YES;
1737	}
1738	else if (expecting == Expecting(ALPHA))
1739	{
1740	strcpy (comptreeParams.comptFileName1, tkn);
1741	MrBayesPrint ("%s Setting comparetree filename 1 to %s\n", spacer, comptreeParams.comptFileName1);
1742	expecting = Expecting(PARAMETER) \| Expecting(SEMICOLON);
1743	}
1744	else
1745	return (ERROR);
1746	}
1747	/* set Filename (comptreeParams.comptFileName2) *************************************************/
1748	else if (!strcmp(parmName, "Filename2"))
1749	{
1750	if (expecting == Expecting(EQUALSIGN))
1751	{
1752	expecting = Expecting(ALPHA);
1753	readWord = YES;
1754	}
1755	else if (expecting == Expecting(ALPHA))
1756	{
1757	strcpy (comptreeParams.comptFileName2, tkn);
1758	MrBayesPrint ("%s Setting comparetree filename 2 to %s\n", spacer, comptreeParams.comptFileName2);
1759	expecting = Expecting(PARAMETER) \| Expecting(SEMICOLON);
1760	}
1761	else
1762	return (ERROR);
1763	}
1764	/* set Filename (comptreeParams.comptOutfile) ***************************************************/
1765	else if (!strcmp(parmName, "Outputname"))
1766	{
1767	if (expecting == Expecting(EQUALSIGN))
1768	{
1769	expecting = Expecting(ALPHA);
1770	readWord = YES;
1771	}
1772	else if (expecting == Expecting(ALPHA))
1773	{
1774	strcpy (comptreeParams.comptOutfile, tkn);
1775	MrBayesPrint ("%s Setting comparetree output file to %s\n", spacer, comptreeParams.comptOutfile);
1776	expecting = Expecting(PARAMETER) \| Expecting(SEMICOLON);
1777	}
1778	else
1779	return (ERROR);
1780	}
1781	/* set Relburnin (chainParams.relativeBurnin) ********************************************************/
1782	else if (!strcmp(parmName, "Relburnin"))
1783	{
1784	if (expecting == Expecting(EQUALSIGN))
1785	expecting = Expecting(ALPHA);
1786	else if (expecting == Expecting(ALPHA))
1787	{
1788	if (IsArgValid(tkn, tempStr) == NO_ERROR)
1789	{
1790	if (!strcmp(tempStr, "Yes"))
1791	chainParams.relativeBurnin = YES;
1792	else
1793	chainParams.relativeBurnin = NO;
1794	}
1795	else
1796	{
1797	MrBayesPrint ("%s Invalid argument for Relburnin\n", spacer);
1798	return (ERROR);
1799	}
1800	if (chainParams.relativeBurnin == YES)
1801	MrBayesPrint ("%s Using relative burnin (a fraction of samples discarded).\n", spacer);
1802	else
1803	MrBayesPrint ("%s Using absolute burnin (a fixed number of samples discarded).\n", spacer);
1804	expecting = Expecting(PARAMETER) \| Expecting(SEMICOLON);
1805	}
1806	else
1807	{
1808	//free (tempStr);
1809	return (ERROR);
1810	}
1811	}
1812	/* set Burnin (chainParams.chainBurnIn) *******************************************************/
1813	else if (!strcmp(parmName, "Burnin"))
1814	{
1815	if (expecting == Expecting(EQUALSIGN))
1816	expecting = Expecting(NUMBER);
1817	else if (expecting == Expecting(NUMBER))
1818	{
1819	sscanf (tkn, "%d", &tempI);
1820	chainParams.chainBurnIn = tempI;
1821	MrBayesPrint ("%s Setting burnin to %ld\n", spacer, chainParams.chainBurnIn);
1822	expecting = Expecting(PARAMETER) \| Expecting(SEMICOLON);
1823	}
1824	else
1825	return (ERROR);
1826	}
1827	/* set Burninfrac (chainParams.burninFraction) ************************************************************/
1828	else if (!strcmp(parmName, "Burninfrac"))
1829	{
1830	if (expecting == Expecting(EQUALSIGN))
1831	expecting = Expecting(NUMBER);
1832	else if (expecting == Expecting(NUMBER))
1833	{
1834	sscanf (tkn, "%lf", &tempD);
1835	if (tempD < 0.01)
1836	{
1837	MrBayesPrint ("%s Burnin fraction too low (< 0.01)\n", spacer);
1838	return (ERROR);
1839	}
1840	if (tempD > 0.50)
1841	{
1842	MrBayesPrint ("%s Burnin fraction too high (> 0.50)\n", spacer);
1843	return (ERROR);
1844	}
1845	chainParams.burninFraction = tempD;
1846	MrBayesPrint ("%s Setting burnin fraction to %.2f\n", spacer, chainParams.burninFraction);
1847	expecting = Expecting(PARAMETER) \| Expecting(SEMICOLON);
1848	}
1849	else
1850	{
1851	return (ERROR);
1852	}
1853	}
1854	/* set Minpartfreq (comptreeParams.minPartFreq) *******************************************************/
1855	else if (!strcmp(parmName, "Minpartfreq"))
1856	{
1857	if (expecting == Expecting(EQUALSIGN))
1858	expecting = Expecting(NUMBER);
1859	else if (expecting == Expecting(NUMBER))
1860	{
1861	sscanf (tkn, "%lf", &tempD);
1862	comptreeParams.minPartFreq = tempD;
1863	MrBayesPrint ("%s Including partitions with probability greater than or equal to %lf in summary statistics\n", spacer, comptreeParams.minPartFreq);
1864	expecting = Expecting(PARAMETER) \| Expecting(SEMICOLON);
1865	}
1866	else
1867	return (ERROR);
1868	}
1869	else
1870	return (ERROR);
1871	}
1872
1873	return (NO_ERROR);
1874
1875	}
1876
1877
1878
1879	#if defined (PRINT_RATEMULTIPLIERS_CPP)
1880	int DELETE_ME_count_taxa(PolyNode *p)
1881	{
1882	int sum=0;
1883
1884	if(p->left==NULL){
1885	if( p->depth > 0.1 )
1886	return 1;
1887	else
1888	return 0;
1889	}
1890
1891	p=p->left;
1892	while(p != NULL){
1893	sum+=DELETE_ME_count_taxa(p);
1894	p=p->sib;
1895	}
1896	return sum;
1897	}
1898
1899
1900
1901
1902	void DELETE_ME_dump_depth(PolyNode *p)
1903	{
1904
1905	/print depth of two taxa clade/
1906	/*if( p->left != NULL && p->left->left == NULL && p->left->sib != NULL && p->left->sib->left == NULL ){
1907	fprintf(rateMultfp,"%f\n",p->depth);
1908	}
1909	*/
1910	/*
1911	if( p->left != NULL && p->left->left == NULL && p->left->sib != NULL && p->left->sib->left == NULL ){
1912	if( p->left->depth > 0.1 && p->left->sib->depth > 0.1 )
1913	fprintf(rateMultfp,"%f\n",p->depth);
1914	}
1915	*/
1916	/print depth of three taxa clade/
1917	if( ((p->left != NULL && p->left->left == NULL) && p->left->sib != NULL && p->left->sib->left != NULL && p->left->sib->left->left == NULL && p->left->sib->left->sib->left == NULL) \|\|
1918	(p->left != NULL && p->left->left != NULL && p->left->left->left == NULL && p->left->left->sib->left == NULL && (p->left->sib->left == NULL)) ){
1919	if( DELETE_ME_count_taxa(p)==2 )
1920	fprintf(rateMultfp,"%f\n",p->depth);
1921	}
1922
1923	p=p->left;
1924	while(p != NULL){
1925	DELETE_ME_dump_depth(p);
1926	p=p->sib;
1927	}
1928	}
1929	#endif
1930
1931
1932
1933
1934
1935	int DoSumt (void)
1936
1937	{
1938
1939	int i, j=0, k, n, len, longestName, treeNo, numTreePartsToPrint,
1940	maxWidthID, maxWidthNumberPartitions, maxNumTaxa, tableWidth=0, unreliable, oneUnreliable,
1941	longestHeader;
1942	MrBFlt f, var_s, sum_s, stddev_s=0.0, sumsq_s, sumStdDev=0.0, maxStdDev=0.0, sumPSRF=0.0,
1943	maxPSRF=0.0, avgStdDev=0.0, avgPSRF=0.0, min_s=0.0, max_s=0.0, numPSRFSamples=0, min;
1944	PartCtr *x;
1945	char *s=NULL, tempName[120], fileName[120], treeName[100], divString[100];
1946	char tempStr=NULL; /not static because error ext is handeled*/
1947	int tempStrLength;
1948	FILE *fp=NULL;
1949	PartCtr *treeParts=NULL,tmp;
1950	SumtFileInfo sumtFileInfo;
1951	Stat theStats;
1952	SafeLong *mask;
1953
1954	#define SCREEN_WIDTH 80
1955
1956	# if defined (MPI_ENABLED)
1957	if (proc_id == 0)
1958	{
1959	# endif
1960
1961
1962	/* Ensure that we read trees with sumt code and not user tree code */
1963	inSumtCommand = YES;
1964
1965	/* set file pointers to NULL */
1966	fp = fpParts = fpTstat = fpVstat = fpCon = fpTrees = NULL;
1967
1968	strcpy(treeName,"tree"); //in case if paramiter is not specified in a .t file
1969
1970	/* Check if there is anything to do */
1971	if (sumtParams.table == NO && sumtParams.summary == NO && sumtParams.showConsensus == NO)
1972	{
1973	MrBayesPrint ("%s Nothing to do, all output parameters (Table, Summary, Consensus) set to 'NO'\n", spacer);
1974	goto errorExit;
1975	}
1976
1977	SafeStrcpy(&tempStr, " ");
1978
1979	/* Initialize sumtParams struct */
1980	sumtParams.numTaxa = 0;
1981	sumtParams.taxaNames = NULL;
1982	sumtParams.SafeLongsNeeded = 0;
1983	sumtParams.tree = AllocatePolyTree (numTaxa);
1984	sumtParams.nESets = 0;
1985	sumtParams.nBSets = 0;
1986	sumtParams.eSetName = NULL;
1987	sumtParams.bSetName = NULL;
1988	sumtParams.popSizeSet = NO;
1989	sumtParams.popSizeSetName = NULL;
1990	AllocatePolyTreePartitions (sumtParams.tree);
1991	sumtParams.numFileTrees = (int ) SafeCalloc (2sumtParams.numRuns+2*numTaxa, sizeof(int));
1992	sumtParams.numFileTreesSampled = sumtParams.numFileTrees + sumtParams.numRuns;
1993	sumtParams.order = sumtParams.numFileTrees + 2*sumtParams.numRuns;
1994	sumtParams.absentTaxa = sumtParams.numFileTrees + 2*sumtParams.numRuns + numTaxa;
1995	if (!sumtParams.numFileTrees)
1996	{
1997	MrBayesPrint ("%s Problems allocating sumtParams.numFileTrees in DoSumt()\n", spacer);
1998	goto errorExit;
1999	}
2000	else
2001	memAllocs[ALLOC_SUMTPARAMS] = YES;
2002
2003	for (treeNo = 0; treeNo < sumtParams.numTrees; treeNo++)
2004	{
2005	/* initialize across-file tree and partition counters */
2006	sumtParams.numTreesSampled = sumtParams.numTreesEncountered = 0;
2007	numUniqueSplitsFound = numUniqueTreesFound = 0;
2008
2009	/* initialize oneUnreliable && unreliable */
2010	oneUnreliable = unreliable = NO;
2011
2012	/* initialize summary statistics */
2013	sumStdDev = 0.0;
2014	sumPSRF = 0.0;
2015	numPSRFSamples = 0;
2016	maxStdDev = 0.0;
2017	maxPSRF = 0.0;
2018
2019	/* tell user we are ready to go */
2020	if (sumtParams.numTrees > 1)
2021	sprintf (fileName,"%s.tree%d", sumtParams.sumtFileName, treeNo+1);
2022	else
2023	strcpy (fileName, sumtParams.sumtFileName);
2024
2025	if (sumtParams.numRuns == 1)
2026	MrBayesPrint ("%s Summarizing trees in file \"%s.t\"\n", spacer, fileName);
2027	else if (sumtParams.numRuns == 2)
2028	MrBayesPrint ("%s Summarizing trees in files \"%s.run1.t\" and \"%s.run2.t\"\n", spacer, fileName, fileName);
2029	else if (sumtParams.numRuns > 2)
2030	MrBayesPrint ("%s Summarizing trees in files \"%s.run1.t\", \"%s.run2.t\",...,\"%s.run%d.t\"\n", spacer, fileName, fileName,fileName,sumtParams.numRuns);
2031
2032	if (chainParams.relativeBurnin == YES)
2033	MrBayesPrint ("%s Using relative burnin ('relburnin=yes'), discarding the first %.0f %% of sampled trees\n",
2034	spacer, chainParams.burninFraction*100.0, chainParams.burninFraction);
2035	else
2036	MrBayesPrint ("%s Using absolute burnin ('relburnin=no'), discarding the first %d sampled trees\n",
2037	spacer, chainParams.chainBurnIn, chainParams.chainBurnIn);
2038
2039	MrBayesPrint ("%s Writing statistics to files %s.<parts\|tstat\|vstat\|trprobs\|con>\n", spacer, sumtParams.sumtOutfile);
2040
2041	for (sumtParams.runId=0; sumtParams.runId < sumtParams.numRuns; sumtParams.runId++)
2042	{
2043	/* initialize tree counters */
2044	sumtParams.numFileTrees[sumtParams.runId] = sumtParams.numFileTreesSampled[sumtParams.runId] = 0;
2045
2046	/* open binary file */
2047	if (sumtParams.numRuns == 1)
2048	sprintf (tempName, "%s.t", fileName);
2049	else
2050	sprintf (tempName, "%s.run%d.t", fileName, sumtParams.runId+1);
2051	strcpy(sumtParams.curFileName, tempName);
2052
2053	/* tell user we are examining files if for the first run */
2054	if (sumtParams.runId == 0)
2055	{
2056	if (sumtParams.numRuns > 1 && sumtParams.numTrees > 1)
2057	MrBayesPrint ("%s Examining first file for tree %d ...\n", spacer, treeNo);
2058	else if (sumtParams.numRuns > 1 && sumtParams.numTrees == 1)
2059	MrBayesPrint ("%s Examining first file ...\n", spacer);
2060	else if (sumtParams.numRuns == 1 && sumtParams.numTrees > 1)
2061	MrBayesPrint ("%s Examining file for tree %d ...\n", spacer, treeNo);
2062	else
2063	MrBayesPrint ("%s Examining file ...\n", spacer);
2064	}
2065
2066	/* examine file */
2067	if (ExamineSumtFile(tempName, &sumtFileInfo, treeName, &sumtParams.brlensDef) == ERROR)
2068	goto errorExit;
2069
2070	/* capture values */
2071	if (sumtParams.runId == 0)
2072
2073	/* catch lack of sampled trees */
2074	if (chainParams.relativeBurnin == NO && chainParams.chainBurnIn > sumtFileInfo.numTreesInLastBlock)
2075	{
2076	MrBayesPrint ("%s No trees are sampled as the burnin exceeds the number of trees in last block\n", spacer);
2077	MrBayesPrint ("%s Try setting burnin to a number less than %d\n", spacer, sumtFileInfo.numTreesInLastBlock);
2078	goto errorExit;
2079	}
2080
2081	/* tell the user that everything is fine */
2082	if (sumtParams.runId == 0)
2083	{
2084	if (sumtFileInfo.numTreeBlocks == 1)
2085	MrBayesPrint ("%s Found one tree block in file \"%s\" with %d trees in last block\n",
2086	spacer, tempName, sumtFileInfo.numTreesInLastBlock);
2087	else
2088	{
2089	MrBayesPrint ("%s Found %d tree blocks in file \"%s\" with %d trees in last block\n",
2090	spacer, sumtFileInfo.numTreeBlocks, tempName, sumtFileInfo.numTreesInLastBlock);
2091	MrBayesPrint ("%s Only the %d trees in last tree block will be summarized\n", spacer, sumtFileInfo.numTreesInLastBlock);
2092	}
2093	sumtParams.numTreesInLastBlock = sumtFileInfo.numTreesInLastBlock;
2094	if (sumtParams.numRuns > 1)
2095	MrBayesPrint ("%s Expecting the same number of trees in the last tree block of all files\n", spacer);
2096	if (chainParams.relativeBurnin == NO)
2097	sumtParams.burnin = chainParams.chainBurnIn;
2098	else
2099	sumtParams.burnin = (int) (sumtFileInfo.numTreesInLastBlock * chainParams.burninFraction);
2100	}
2101	else
2102	{
2103	if (sumtFileInfo.numTreesInLastBlock != sumtParams.numFileTrees[0])
2104	{
2105	MrBayesPrint ("%s Found %d trees in first file but %d trees in file \"%s\"\n", spacer,
2106	sumtParams.numFileTrees[0],
2107	sumtFileInfo.numTreesInLastBlock,
2108	tempName);
2109	goto errorExit;
2110	}
2111	}
2112
2113	/* Now we read the file for real. First, allocate a string for reading the file... */
2114	if (sumtParams.runId == 0 && treeNo == 0)
2115	{
2116	s = (char )SafeMalloc((size_t) (sumtFileInfo.longestLineLength sizeof(char)));
2117	if (!s)
2118	{
2119	MrBayesPrint ("%s Problem allocating string for reading sumt file\n", spacer);
2120	goto errorExit;
2121	}
2122	}
2123	else
2124	{
2125	free (s);
2126	s = (char ) SafeMalloc (sizeof (char) sumtFileInfo.longestLineLength);
2127	if (!s)
2128	{
2129	MrBayesPrint ("%s Problem reallocating string for reading sumt file\n", spacer);
2130	goto errorExit;
2131	}
2132	}
2133
2134	/* ... open the file ... */
2135	if ((fp = OpenTextFileR(tempName)) == NULL)
2136	goto errorExit;
2137
2138	/* ...and fast forward to beginning of last tree block. */
2139	for (i=0; i<sumtFileInfo.lastTreeBlockBegin+1; i++)
2140	{
2141	if( fgets (s, sumtFileInfo.longestLineLength-2, fp) == NULL )
2142	{
2143	printf("Error in function: %s at line: %d in file: %s", __FUNCTION__, __LINE__, __FILE__);
2144	}
2145	}
2146
2147
2148	#if defined (PRINT_RATEMULTIPLIERS_CPP)
2149	sprintf (tempName, "%s.ratemult", chainParams.chainFileName);
2150	if ( (rateMultfp=OpenNewMBPrintFile (tempName)) == NULL )
2151	{
2152	printf("Error oppening file: %s to write", tempName);
2153	goto errorExit;
2154	}
2155	fprintf(rateMultfp,"rateMult_CPP\n");
2156	#endif
2157
2158
2159	/* Set up cheap status bar. */
2160	if (sumtParams.runId == 0)
2161	{
2162	if (sumtParams.numTrees > 1)
2163	MrBayesPrint ("\n%s Tree reading status for tree %d:\n\n", spacer, treeNo+1);
2164	else
2165	MrBayesPrint ("\n%s Tree reading status:\n\n", spacer);
2166	MrBayesPrint ("%s 0 10 20 30 40 50 60 70 80 90 100\n", spacer);
2167	MrBayesPrint ("%s v-------v-------v-------v-------v-------v-------v-------v-------v-------v-------v\n", spacer);
2168	MrBayesPrint ("%s *", spacer);
2169	numAsterices = 0;
2170	}
2171
2172	/* ...and parse file, tree-by-tree. We are only parsing lines between the "begin trees" and "end" statements.
2173	We don't actually get those lines, however, but rather the lines between those statements. */
2174	expecting = Expecting(COMMAND);
2175	/* We skip the begin trees statement so we need to set up some variables here */
2176	inTreesBlock = YES;
2177	ResetTranslateTable();
2178	for (i=0; i<sumtFileInfo.lastTreeBlockEnd - sumtFileInfo.lastTreeBlockBegin - 1; i++)
2179	{
2180	if( fgets (s, sumtFileInfo.longestLineLength-2, fp) == NULL )
2181	{
2182	printf("Error in function: %s at line: %d in file: %s", __FUNCTION__, __LINE__, __FILE__);
2183	}
2184	/MrBayesPrint ("%s", s);/
2185	if (ParseCommand (s) == ERROR)
2186	goto errorExit;
2187	}
2188	inTreesBlock = NO;
2189	ResetTranslateTable();
2190
2191	/* Finish cheap status bar. */
2192	if (sumtParams.runId == sumtParams.numRuns - 1)
2193	{
2194	if (numAsterices < 80)
2195	{
2196	for (i=0; i<80 - numAsterices; i++)
2197	MrBayesPrint ("*");
2198	}
2199	MrBayesPrint ("\n\n");
2200	}
2201
2202	/* print out information on absent and pruned taxa */
2203	if (sumtParams.runId == sumtParams.numRuns - 1 && treeNo == 0)
2204	PrintSumtTaxaInfo ();
2205
2206	/* tell user how many trees were successfully read */
2207	if (sumtParams.numRuns == 1)
2208	MrBayesPrint ("%s Read %d trees from last tree block (sampling %d of them)\n", spacer,
2209	sumtParams.numTreesEncountered, sumtParams.numTreesSampled);
2210	else if (sumtParams.numRuns > 1)
2211	{
2212	if (sumtParams.runId != 0 && sumtParams.numFileTreesSampled[sumtParams.runId]!=sumtParams.numFileTreesSampled[0])
2213	{
2214	if (sumtParams.runId == sumtParams.numRuns - 1)
2215	MrBayesPrint ("\n\n");
2216	MrBayesPrint ("%s Found %d post-burnin trees in the first file but %d post-burnin trees in file %d\n",
2217	spacer,
2218	sumtParams.numFileTreesSampled[0],
2219	sumtParams.numFileTreesSampled[sumtParams.runId],
2220	sumtParams.runId+1);
2221	goto errorExit;
2222	}
2223	if (sumtParams.runId == sumtParams.numRuns - 1)
2224	{
2225	MrBayesPrint ("%s Read a total of %d trees in %d files (sampling %d of them)\n", spacer,
2226	sumtParams.numTreesEncountered,
2227	sumtParams.numRuns, sumtParams.numTreesSampled);
2228	MrBayesPrint ("%s (Each file contained %d trees of which %d were sampled)\n", spacer,
2229	sumtParams.numFileTrees[0],
2230	sumtParams.numFileTreesSampled[0]);
2231	}
2232	}
2233
2234	/* Check that at least one tree was read in. */
2235	if (sumtParams.numTreesSampled <= 0)
2236	{
2237	MrBayesPrint ("%s No trees read in\n", spacer);
2238	goto errorExit;
2239	}
2240
2241	SafeFclose (&fp);
2242	} /* next run for this tree */
2243
2244	/* Extract partition counter pointers */
2245	treeParts = (PartCtr *) SafeCalloc ((size_t)(numUniqueSplitsFound), sizeof(PartCtr ));
2246	i = 0;
2247	PartCtrUppass(partCtrRoot, treeParts, &i);
2248
2249	min = (sumtParams.minPartFreq * (sumtParams.numTreesSampled/sumtParams.numRuns));
2250	numTreePartsToPrint=numUniqueSplitsFound;
2251	for (i=0; i<numTreePartsToPrint;)
2252	{
2253	for (j=0; j<sumtParams.numRuns;j++)
2254	{
2255	if(treeParts[i]->count[j]>=min)
2256	break;
2257	}
2258	if( j==sumtParams.numRuns )
2259	{
2260	numTreePartsToPrint--;
2261	tmp=treeParts[numTreePartsToPrint];
2262	treeParts[numTreePartsToPrint]=treeParts[i];
2263	treeParts[i]=tmp;
2264	}
2265	else
2266	{
2267	i++;
2268	}
2269	}
2270
2271
2272	/* Sort taxon partitions (clades, splits) ... */
2273	SortPartCtr (treeParts, 0, numTreePartsToPrint-1);
2274
2275	/* Sort root and tips among those splits always present */
2276	SortTerminalPartCtr (treeParts, numUniqueSplitsFound);
2277
2278	/* open output files for summary information (three files) */
2279	if (OpenSumtFiles (treeNo) == ERROR)
2280	goto errorExit;
2281
2282	/* Print partitions to screen. */
2283	if (treeNo == 0)
2284	{
2285	longestName = 0;
2286	for (k=0; k<numTaxa; k++)
2287	{
2288	if (taxaInfo[k].isDeleted == NO && sumtParams.absentTaxa[k] == NO)
2289	continue;
2290	len = (int) strlen (taxaNames[k]);
2291	if (len > longestName)
2292	longestName = len;
2293	}
2294	if (sumtParams.table == YES)
2295	{
2296	MrBayesPrint (" \n");
2297	MrBayesPrint ("%s General explanation: \n", spacer);
2298	MrBayesPrint (" \n");
2299	MrBayesPrint ("%s In an unrooted tree, a taxon bipartition (split) is specified by removing a \n", spacer);
2300	MrBayesPrint ("%s branch, thereby dividing the species into those to the left and those to the \n", spacer);
2301	MrBayesPrint ("%s right of the branch. Here, taxa to one side of the removed branch are denoted \n", spacer);
2302	MrBayesPrint ("%s '.' and those to the other side are denoted '*'. Specifically, the '.' symbol \n", spacer);
2303	MrBayesPrint ("%s is used for the taxa on the same side as the outgroup. \n", spacer);
2304	MrBayesPrint (" \n");
2305	MrBayesPrint ("%s In a rooted or clock tree, the tree is rooted using the model and not by \n", spacer);
2306	MrBayesPrint ("%s reference to an outgroup. Each bipartition therefore corresponds to a clade, \n", spacer);
2307	MrBayesPrint ("%s that is, a group that includes all the descendants of a particular branch in \n", spacer);
2308	MrBayesPrint ("%s the tree. Taxa that are included in each clade are denoted using '*', and \n", spacer);
2309	MrBayesPrint ("%s taxa that are not included are denoted using the '.' symbol. \n", spacer);
2310	MrBayesPrint (" \n");
2311	MrBayesPrint ("%s The output first includes a key to all the bipartitions with frequency larger \n", spacer);
2312	MrBayesPrint ("%s or equual to (Minpartfreq) in at least one run. Minpartfreq is a paramiter to \n", spacer);
2313	MrBayesPrint ("%s sumt command and currently it is set to %1.2lf. This is followed by a table \n", spacer, sumtParams.minPartFreq);
2314	MrBayesPrint ("%s with statistics for the informative bipartitions (those including at least \n", spacer);
2315	MrBayesPrint ("%s two taxa), sorted from highest to lowest probability. For each bipartition, \n", spacer);
2316	MrBayesPrint ("%s the table gives the number of times the partition or split was observed in all\n", spacer);
2317	MrBayesPrint ("%s runs (#obs) and the posterior probability of the bipartition (Probab.), which \n", spacer);
2318	MrBayesPrint ("%s is the same as the split frequency. If several runs are summarized, this is \n", spacer);
2319	MrBayesPrint ("%s followed by the minimum split frequency (Min(s)), the maximum frequency \n", spacer);
2320	MrBayesPrint ("%s (Max(s)), and the standard deviation of frequencies (Stddev(s)) across runs. \n", spacer);
2321	MrBayesPrint ("%s The latter value should approach 0 for all bipartitions as MCMC runs converge.\n", spacer);
2322	MrBayesPrint (" \n");
2323	MrBayesPrint ("%s This is followed by a table summarizing branch lengths, node heights (if a \n", spacer);
2324	MrBayesPrint ("%s clock model was used) and relaxed clock parameters (if a relaxed clock model \n", spacer);
2325	MrBayesPrint ("%s was used). The mean, variance, and 95 %% credible interval are given for each \n", spacer);
2326	MrBayesPrint ("%s of these parameters. If several runs are summarized, the potential scale \n", spacer);
2327	MrBayesPrint ("%s reduction factor (PSRF) is also given; it should approach 1 as runs converge. \n", spacer);
2328	MrBayesPrint ("%s Node heights will take calibration points into account, if such points were \n", spacer);
2329	MrBayesPrint ("%s used in the analysis. \n", spacer);
2330	MrBayesPrint ("%s \n", spacer);
2331	MrBayesPrint ("%s Note that Stddev may be unreliable if the partition is not present in all \n", spacer);
2332	MrBayesPrint ("%s runs (the last column indicates the number of runs that sampled the partition \n", spacer);
2333	MrBayesPrint ("%s if more than one run is summarized). The PSRF is not calculated at all if \n", spacer);
2334	MrBayesPrint ("%s the partition is not present in all runs.The PSRF is also sensitive to small \n", spacer);
2335	MrBayesPrint ("%s sample sizes and it should only be considered a rough guide to convergence \n", spacer);
2336	MrBayesPrint ("%s since some of the assumptions allowing one to interpret it as a true potential\n", spacer);
2337	MrBayesPrint ("%s scale reduction factor are violated in MrBayes. \n", spacer);
2338	MrBayesPrint ("%s \n", spacer);
2339	MrBayesPrint ("%s List of taxa in bipartitions: \n", spacer);
2340	MrBayesPrint (" \n");
2341	j = 1;
2342	for (k=0; k<numTaxa; k++)
2343	{
2344	if (taxaInfo[k].isDeleted == NO && sumtParams.absentTaxa[k] == NO)
2345	{
2346	MrBayesPrint ("%s %4d -- %s\n", spacer, j++, taxaNames[k]);
2347	}
2348	}
2349	}
2350	}
2351
2352	if (sumtParams.numTrees > 1 && (sumtParams.table == YES \|\| sumtParams.summary == YES))
2353	{
2354	MrBayesPrint ("\n\n");
2355	MrBayesPrint ("%s Results for tree number %d\n", spacer, treeNo+1);
2356	MrBayesPrint ("%s ==========================\n\n", spacer);
2357	}
2358
2359	/* First print key to taxon bipartitions */
2360	if (sumtParams.table == YES)
2361	{
2362	MrBayesPrint ("\n");
2363	if (sumtParams.numTrees == 1)
2364	MrBayesPrint ("%s Key to taxon bipartitions (saved to file \"%s.parts\"):\n\n", spacer, sumtParams.sumtOutfile);
2365	else
2366	MrBayesPrint ("%s Key to taxon bipartitions (saved to file \"%s.tree%d.parts\"):\n\n", spacer, sumtParams.sumtOutfile, treeNo+1 );
2367	}
2368
2369	/* calculate a couple of numbers that are handy to have */
2370	/*numTreePartsToPrint = 0;
2371	for (i=0; i<numUniqueSplitsFound; i++)
2372	{
2373	if ((MrBFlt)treeParts[i]->totCount/(MrBFlt)sumtParams.numTreesSampled >= sumtParams.minPartFreq)
2374	numTreePartsToPrint++;
2375	}
2376	*/
2377	maxWidthID = (int) (log10 (numTreePartsToPrint)) + 1;
2378	if (maxWidthID < 2)
2379	maxWidthID = 2;
2380	maxNumTaxa = SCREEN_WIDTH - 9;
2381
2382	for (j=0; j<sumtParams.numTaxa; j+=maxNumTaxa)
2383	{
2384	/* print header to screen and to parts file simultaneously */
2385
2386	/* first print header to screen */
2387	MrBayesPrint ("%s ", spacer);
2388	if (j == 0)
2389	MrBayesPrint ("%*s -- Partition\n", maxWidthID, "ID");
2390	else
2391	MrBayesPrint ("%*s -- Partition (continued)\n", maxWidthID, "ID");
2392	tableWidth = maxWidthID + 4 + sumtParams.numTaxa;
2393	if (tableWidth > SCREEN_WIDTH)
2394	tableWidth = SCREEN_WIDTH;
2395	MrBayesPrint ("%s ", spacer);
2396	for (i=0; i<tableWidth; i++)
2397	MrBayesPrint ("-");
2398	MrBayesPrint ("\n");
2399
2400	/* now print header to file */
2401	MrBayesPrintf (fpParts, "ID\tPartition\n");
2402
2403	/* now, show partitions that were found on screen; print to .parts file simultaneously */
2404	mask = SafeCalloc (sumtParams.SafeLongsNeeded, sizeof(SafeLong));
2405	for (i=0; i<sumtParams.numTaxa; i++)
2406	SetBit (i, mask);
2407	for (i=0; i<numTreePartsToPrint; i++)
2408	{
2409	x = treeParts[i];
2410	if (IsBitSet(localOutGroup, x->partition) == YES && sumtParams.isRooted == NO)
2411	FlipBits(x->partition, sumtParams.SafeLongsNeeded, mask);
2412
2413	if ((NumBits(x->partition, sumtParams.SafeLongsNeeded) == numLocalTaxa \|\| NumBits(x->partition, sumtParams.SafeLongsNeeded) == 0) && sumtParams.isClock == NO)
2414	continue;
2415
2416	MrBayesPrint ("%s %*d -- ", spacer, maxWidthID, i);
2417	if (sumtParams.numTaxa <= maxNumTaxa)
2418	ShowParts (stdout, x->partition, sumtParams.numTaxa);
2419	else
2420	{
2421	if (sumtParams.numTaxa - j > maxNumTaxa)
2422	ShowSomeParts (stdout, x->partition, j, maxNumTaxa);
2423	else
2424	ShowSomeParts (stdout, x->partition, j, sumtParams.numTaxa - j);
2425	}
2426	fflush(stdout);
2427	MrBayesPrint ("\n");
2428
2429	MrBayesPrintf (fpParts, "%d\t", i);
2430	ShowParts (fpParts, x->partition, sumtParams.numTaxa);
2431	MrBayesPrintf (fpParts, "\n");
2432	}
2433	free (mask);
2434
2435	/* finish screen table */
2436	if (sumtParams.table == YES)
2437	{
2438	MrBayesPrint ("%s ", spacer);
2439	for (i=0; i<tableWidth; i++)
2440	{
2441	MrBayesPrint ("-");
2442	}
2443	MrBayesPrint ("\n");
2444	if (oneUnreliable == YES)
2445	{
2446	MrBayesPrint ("%s * The partition was not found in all runs so the values are unreliable\n\n", spacer);
2447	}
2448	else
2449	{
2450	MrBayesPrint ("\n");
2451	}
2452	}
2453	}
2454
2455	/* Second, print statitistics for taxon bipartitions */
2456	if (sumtParams.table == YES)
2457	{
2458	if (sumtParams.isRooted == NO)
2459	MrBayesPrint ("%s Summary statistics for informative taxon bipartitions\n", spacer);
2460	else
2461	MrBayesPrint ("%s Summary statistics for informative taxon bipartitions (clades)\n", spacer);
2462	MrBayesPrint ("%s (saved to file \"%s.tstat\"):\n\n", spacer, sumtParams.sumtOutfile);
2463	}
2464
2465	/* calculate a couple of numbers that are handy to have */
2466	/*numTreePartsToPrint = 0;
2467	for (i=0; i<numUniqueSplitsFound; i++)
2468	{
2469	if ((MrBFlt)treeParts[i]->totCount/(MrBFlt)sumtParams.numTreesSampled >= sumtParams.minPartFreq)
2470	numTreePartsToPrint++;
2471	}
2472	*/
2473	maxWidthID = (int) (log10 (numTreePartsToPrint)) + 1;
2474	if (maxWidthID < 2)
2475	maxWidthID = 2;
2476	maxWidthNumberPartitions = (int) (log10 (treeParts[0]->totCount)) + 1;
2477	if (maxWidthNumberPartitions < 4)
2478	maxWidthNumberPartitions = 4;
2479
2480	/* print header to screen and to parts file simultaneously */
2481	if (sumtParams.table == YES)
2482	{
2483	/* first print header to screen */
2484	MrBayesPrint ("%s ", spacer);
2485	MrBayesPrint ("%*s ", maxWidthID, "ID");
2486	tableWidth = maxWidthID + 3;
2487	MrBayesPrint ("#obs");
2488	tableWidth += 4;
2489	for (i=4; i<maxWidthNumberPartitions; i++)
2490	{
2491	MrBayesPrint (" ");
2492	tableWidth++;
2493	}
2494	MrBayesPrint (" Probab.");
2495	tableWidth += 11;
2496	if (sumtParams.numRuns > 1)
2497	{
2498	MrBayesPrint (" Sd(s)+ ");
2499	MrBayesPrint (" Min(s) Max(s) ");
2500	tableWidth += 36;
2501	MrBayesPrint (" Nruns ");
2502	tableWidth += 8;
2503	}
2504	MrBayesPrint ("\n%s ", spacer);
2505	for (i=0; i<tableWidth; i++)
2506	{
2507	MrBayesPrint ("-");
2508	}
2509	MrBayesPrint ("\n");
2510
2511	/* now print header to file */
2512	if (sumtParams.numRuns > 1)
2513	MrBayesPrintf (fpTstat, "ID\t#obs\tProbability(=s)\tStddev(s)\tMin(s)\tMax(s)\tNruns\n");
2514	else
2515	MrBayesPrintf (fpTstat, "ID\t#obs\tProbability(=s)\n");
2516	}
2517
2518	/* now, show informative partitions that were found on screen; print to .tstat file simultaneously */
2519	for (i=0; i<numTreePartsToPrint; i++)
2520	{
2521	x = treeParts[i];
2522
2523	/* skip uninformative partitions */
2524	if (NumBits(x->partition, sumtParams.SafeLongsNeeded) <= 1 \|\| NumBits(x->partition, sumtParams.SafeLongsNeeded) == sumtParams.numTaxa)
2525	continue;
2526	if (NumBits(x->partition, sumtParams.SafeLongsNeeded) == sumtParams.numTaxa - 1 && sumtParams.isRooted == NO)
2527	continue;
2528
2529	if (sumtParams.table == YES)
2530	{
2531	MrBayesPrint ("%s %*d", spacer, maxWidthID, i);
2532	fflush(stdout);
2533	MrBayesPrintf (fpTstat, "%d\t", i);
2534	}
2535	if (sumtParams.numRuns > 1)
2536	{
2537	sum_s = 0.0;
2538	sumsq_s = 0.0;
2539	min_s = 1.0;
2540	max_s = 0.0;
2541	for (n=j=0; n<sumtParams.numRuns; n++)
2542	{
2543	if (x->count[n] > 0)
2544	j++;
2545	f = (MrBFlt) x->count[n] / (MrBFlt) sumtParams.numFileTreesSampled[n];
2546	sum_s += f;
2547	sumsq_s += f * f;
2548	if (f < min_s)
2549	min_s = f;
2550	if (f > max_s)
2551	max_s = f;
2552	}
2553	var_s = sumsq_s - sum_s * sum_s / (MrBFlt) sumtParams.numRuns;
2554	var_s /= (sumtParams.numRuns - 1);
2555	if (var_s > 0.0)
2556	stddev_s = sqrt (var_s);
2557	else
2558	stddev_s = 0.0;
2559	if (j == sumtParams.numRuns)
2560	unreliable = NO;
2561	else
2562	{
2563	unreliable = YES;
2564	oneUnreliable = YES;
2565	}
2566	}
2567	if (sumtParams.table == YES)
2568	{
2569	f = (MrBFlt) x->totCount / (MrBFlt) sumtParams.numTreesSampled;
2570	MrBayesPrint (" %*d %1.6lf", maxWidthNumberPartitions, x->totCount, f);
2571	MrBayesPrintf (fpTstat, "\t%d\t%s", x->totCount, MbPrintNum(f));
2572	if (sumtParams.numRuns > 1)
2573	{
2574	MrBayesPrint (" %1.6lf %1.6lf %1.6lf", stddev_s, min_s, max_s);
2575	MrBayesPrint (" %3d", j);
2576	MrBayesPrintf (fpTstat, "\t%s", MbPrintNum(stddev_s));
2577	MrBayesPrintf (fpTstat, "\t%s", MbPrintNum(min_s));
2578	MrBayesPrintf (fpTstat, "\t%s", MbPrintNum(max_s));
2579	MrBayesPrintf (fpTstat, "\t%d", j);
2580	}
2581	MrBayesPrintf (fpTstat, "\n");
2582	if (unreliable == YES)
2583	MrBayesPrint (" *\n");
2584	else
2585	MrBayesPrint ("\n");
2586	sumStdDev += stddev_s;
2587	if (stddev_s > maxStdDev)
2588	maxStdDev = stddev_s;
2589	}
2590	}
2591
2592	/* finish screen table */
2593	if (sumtParams.table == YES)
2594	{
2595	MrBayesPrint ("%s ", spacer);
2596	for (i=0; i<tableWidth; i++)
2597	{
2598	MrBayesPrint ("-");
2599	}
2600	MrBayesPrint ("\n");
2601	if (sumtParams.numRuns > 1)
2602	{
2603	MrBayesPrint ("%s + Convergence diagnostic (standard deviation of split frequencies)\n", spacer);
2604	MrBayesPrint ("%s should approach 0.0 as runs converge.\n\n", spacer);
2605	}
2606	if (oneUnreliable == YES)
2607	MrBayesPrint ("%s * The partition was not found in all runs so the values are unreliable\n", spacer);
2608	}
2609
2610	/* Third, print statitistics for branch and node parameters */
2611	if (sumtParams.table == YES)
2612	{
2613	MrBayesPrint ("\n");
2614	MrBayesPrint ("%s Summary statistics for branch and node parameters\n", spacer);
2615	MrBayesPrint ("%s (saved to file \"%s.vstat\"):\n", spacer, sumtParams.sumtOutfile);
2616	}
2617
2618	if (sumtParams.table == YES)
2619	{
2620	/* calculate longest header */
2621	longestHeader = 9; /* length of 'parameter' */
2622	i = (int)(log10(numTreePartsToPrint)) + 3; /* length of partition specifier including [] */
2623	len = i + (int)(strlen(treeName)) + 2; /* length of length{m}[n] or height{m}[n] */
2624	if (len > longestHeader)
2625	longestHeader = len;
2626	for (j=0; j<sumtParams.nBSets; j++)
2627	{
2628	len = (int) strlen(sumtParams.tree->bSetName[j]) + 7 + i;
2629	if (len > longestHeader)
2630	longestHeader = len;
2631	}
2632	for (j=0; j<sumtParams.nESets; j++)
2633	{
2634	len = (int) strlen(sumtParams.tree->eSetName[j]) + 8 + i;
2635	if (len > longestHeader)
2636	longestHeader = len;
2637	}
2638	if (sumtParams.popSizeSet == YES)
2639	{
2640	len = (int) strlen(sumtParams.tree->popSizeSetName) + i;
2641	if (len > longestHeader)
2642	longestHeader = len;
2643	}
2644
2645	/* print the header rows */
2646	MrBayesPrint ("\n");
2647	if (sumtParams.HPD == NO)
2648	MrBayesPrint ("%s %*c 95%% Cred. Interval\n", spacer, longestHeader, ' ');
2649	else
2650	MrBayesPrint ("%s %*c 95%% HPD Interval\n", spacer, longestHeader, ' ');
2651	MrBayesPrint ("%s %*c --------------------\n", spacer, longestHeader, ' ');
2652
2653	MrBayesPrint ("%s Parameter%*c Mean Variance Lower Upper Median", spacer, longestHeader-9, ' ');
2654	tableWidth = 68 + longestHeader - 9;
2655	if (sumtParams.HPD == YES)
2656	MrBayesPrintf (fpVstat, "Parameter\tMean\tVariance\tCredInt_Lower\tCredInt_Upper\tMedian", spacer, longestHeader-9, ' ');
2657	else
2658	MrBayesPrintf (fpVstat, "Parameter\tMean\tVariance\tHPD_Lower\tHPD_Upper\tMedian", spacer, longestHeader-9, ' ');
2659	if (sumtParams.numRuns > 1)
2660	{
2661	MrBayesPrint (" PSRF+ Nruns");
2662	tableWidth += 17;
2663	MrBayesPrintf (fpVstat, "\tPSRF\tNruns");
2664	}
2665	MrBayesPrint ("\n");
2666	MrBayesPrintf (fpVstat, "\n");
2667
2668	MrBayesPrint ("%s ", spacer);
2669	for (j=0; j<tableWidth; j++)
2670	{
2671	MrBayesPrint ("-");
2672	}
2673	MrBayesPrint ("\n");
2674
2675	/* print lengths */
2676	strcpy (divString, treeName+4);
2677	for (i=1; i<numTreePartsToPrint; i++)
2678	{
2679	x = treeParts[i];
2680	tempStrLength=(int)strlen(tempStr);
2681	SafeSprintf (&tempStr,&tempStrLength, "length%s[%d]", divString, i);
2682	len = (int) strlen(tempStr);
2683
2684	GetSummary (x->length, sumtParams.numRuns, x->count, &theStats, sumtParams.HPD);
2685
2686	MrBayesPrint ("%s %-*s ", spacer, longestHeader, tempStr);
2687	MrBayesPrintf (fpVstat, "%s", tempStr);
2688
2689	PrintSumtTableLine(sumtParams.numRuns, x->count, &theStats, &numPSRFSamples, &maxPSRF, &sumPSRF);
2690
2691	}
2692
2693	/* print heights */
2694	if (sumtParams.isClock == YES)
2695	{
2696	strcpy (divString, treeName+4);
2697	for (i=0; i<numTreePartsToPrint; i++)
2698	{
2699	x = treeParts[i];
2700	tempStrLength=(int)strlen(tempStr);
2701	SafeSprintf (&tempStr,&tempStrLength, "height%s[%d]", divString, i);
2702	len = (int) strlen(tempStr);
2703
2704	GetSummary (x->height, sumtParams.numRuns, x->count, &theStats, sumtParams.HPD);
2705
2706	MrBayesPrint ("%s %-*s ", spacer, longestHeader, tempStr);
2707	MrBayesPrintf (fpVstat, "%s", tempStr);
2708
2709	PrintSumtTableLine(sumtParams.numRuns, x->count, &theStats, &numPSRFSamples, &maxPSRF, &sumPSRF);
2710	}
2711	}
2712
2713	/* print ages */
2714	if (sumtParams.isCalibrated == YES)
2715	{
2716	strcpy (divString, treeName+4);
2717	for (i=0; i<numTreePartsToPrint; i++)
2718	{
2719	x = treeParts[i];
2720	tempStrLength=(int)strlen(tempStr);
2721	SafeSprintf (&tempStr,&tempStrLength, "age%s[%d]", divString, i);
2722	len = (int) strlen(tempStr);
2723
2724	GetSummary (x->age, sumtParams.numRuns, x->count, &theStats, sumtParams.HPD);
2725
2726	MrBayesPrint ("%s %-*s ", spacer, longestHeader, tempStr);
2727	MrBayesPrintf (fpVstat, "%s", tempStr);
2728
2729	PrintSumtTableLine(sumtParams.numRuns, x->count, &theStats, &numPSRFSamples, &maxPSRF, &sumPSRF);
2730	}
2731	}
2732
2733	/* print effective branch lengths */
2734	if (sumtParams.isRelaxed == YES)
2735	{
2736	for (i=0; i<sumtParams.nBSets; i++)
2737	{
2738	for (j=1; j<numTreePartsToPrint; j++)
2739	{
2740	x = treeParts[j];
2741	tempStrLength=(int)strlen(tempStr);
2742	SafeSprintf (&tempStr,&tempStrLength, "%s_length[%d]", sumtParams.bSetName[i], j);
2743	len = (int) strlen(tempStr);
2744
2745	GetSummary (x->bLen[i], sumtParams.numRuns, x->count, &theStats, sumtParams.HPD);
2746
2747	MrBayesPrint ("%s %-*s ", spacer, longestHeader, tempStr);
2748	MrBayesPrintf (fpVstat, "%s", tempStr);
2749
2750	PrintSumtTableLine(sumtParams.numRuns, x->count, &theStats, &numPSRFSamples, &maxPSRF, &sumPSRF);
2751	}
2752	for (j=1; j<numTreePartsToPrint; j++)
2753	{
2754	x = treeParts[j];
2755	tempStrLength=(int)strlen(tempStr);
2756	SafeSprintf (&tempStr,&tempStrLength, "%s_rate[%d]", sumtParams.bSetName[i], j);
2757	len = (int) strlen(tempStr);
2758
2759	GetSummary (x->bRate[i], sumtParams.numRuns, x->count, &theStats, sumtParams.HPD);
2760
2761	MrBayesPrint ("%s %-*s ", spacer, longestHeader, tempStr);
2762	MrBayesPrintf (fpVstat, "%s", tempStr);
2763
2764	PrintSumtTableLine(sumtParams.numRuns, x->count, &theStats, &numPSRFSamples, &maxPSRF, &sumPSRF);
2765	}
2766	}
2767	for (i=0; i<sumtParams.nESets; i++)
2768	{
2769	for (j=1; j<numTreePartsToPrint; j++)
2770	{
2771	x = treeParts[j];
2772	tempStrLength=(int)strlen(tempStr);
2773	SafeSprintf (&tempStr,&tempStrLength, "%s_nEvents[%d]", sumtParams.eSetName[i], j);
2774	len = (int) strlen(tempStr);
2775
2776	GetIntSummary (x->nEvents[i], sumtParams.numRuns, x->count, &theStats, sumtParams.HPD);
2777
2778	MrBayesPrint ("%s %-*s ", spacer, longestHeader, tempStr);
2779	MrBayesPrintf (fpVstat, "%s", tempStr);
2780
2781	PrintSumtTableLine(sumtParams.numRuns, x->count, &theStats, &numPSRFSamples, &maxPSRF, &sumPSRF);
2782	}
2783	}
2784	}
2785
2786	/* print population size sets */
2787	if (sumtParams.popSizeSet == YES)
2788	{
2789	for (j=1; j<numTreePartsToPrint; j++)
2790	{
2791	x = treeParts[j];
2792	tempStrLength=(int)strlen(tempStr);
2793	SafeSprintf (&tempStr,&tempStrLength, "%s[%d]", sumtParams.popSizeSetName, j);
2794	len = (int) strlen(tempStr);
2795
2796	GetSummary (x->popSize, sumtParams.numRuns, x->count, &theStats, sumtParams.HPD);
2797
2798	MrBayesPrint ("%s %-*s ", spacer, longestHeader, tempStr);
2799	MrBayesPrintf (fpVstat, "%s", tempStr);
2800
2801	PrintSumtTableLine(sumtParams.numRuns, x->count, &theStats, &numPSRFSamples, &maxPSRF, &sumPSRF);
2802	}
2803	}
2804
2805	/* finish table */
2806	MrBayesPrint ("%s ", spacer);
2807	for (j=0; j<tableWidth; j++)
2808	MrBayesPrint ("-");
2809	MrBayesPrint ("\n");
2810
2811	if (sumtParams.numRuns > 1)
2812	{
2813	MrBayesPrint ("%s + Convergence diagnostic (PSRF = Potential Scale Reduction Factor; Gelman\n", spacer);
2814	MrBayesPrint ("%s and Rubin, 1992) should approach 1.0 as runs converge. NA is reported when\n", spacer);
2815	MrBayesPrint ("%s deviation of parameter values within all runs is 0 or when a parameter\n", spacer);
2816	MrBayesPrint ("%s value (a branch length, for instance) is not sampled in all runs.\n", spacer);
2817	}
2818
2819	if (oneUnreliable == YES)
2820	{
2821	MrBayesPrint ("%s * The partition was not found in all runs so the values are unreliable.\n", spacer);
2822	}
2823	MrBayesPrint ("\n\n");
2824	}
2825
2826	/* Exclude trivial splits when calculating average standard deviation of split frequencies. */
2827	avgStdDev = sumStdDev / (numTreePartsToPrint-sumtParams.numTaxa-1);
2828	avgPSRF = sumPSRF / numPSRFSamples;
2829
2830	if (sumtParams.numRuns > 1 && sumtParams.summary == YES)
2831	{
2832	MrBayesPrint ("%s Summary statistics for partitions with frequency >= %1.2lf in at least one run:\n", spacer, sumtParams.minPartFreq);
2833	MrBayesPrint ("%s Average standard deviation of split frequencies = %1.6lf\n", spacer, avgStdDev);
2834	MrBayesPrint ("%s Maximum standard deviation of split frequencies = %1.6lf\n", spacer, maxStdDev);
2835	}
2836	if (sumtParams.brlensDef == YES && sumtParams.numRuns > 1 && sumtParams.summary == YES)
2837	{
2838	MrBayesPrint ("%s Average PSRF for parameter values ( excluding NA and >10.0 ) = %1.3lf\n", spacer, avgPSRF);
2839	if(maxPSRF == 10 )
2840	MrBayesPrint ("%s Maximum PSRF for parameter values = NA\n", spacer);
2841	else
2842	MrBayesPrint ("%s Maximum PSRF for parameter values = %1.3lf\n", spacer, maxPSRF);
2843	}
2844	MrBayesPrint ("\n");
2845
2846
2847	SortPartCtr (treeParts, 0, numUniqueSplitsFound-1); /* We sort again but this time we sort all partitions instead of just first numTreePartsToPrintNow */
2848	/* make the majority rule consensus tree */
2849	if (sumtParams.showConsensus == YES && ConTree (treeParts, numUniqueSplitsFound) == ERROR)
2850	goto errorExit;
2851
2852	/* get probabilities of individual trees */
2853	if (TreeProb () == ERROR)
2854	goto errorExit;
2855
2856	/* print brlens */
2857	if (sumtParams.printBrlensToFile == YES && PrintBrlensToFile (treeParts, numUniqueSplitsFound, treeNo) == ERROR)
2858	goto errorExit;
2859
2860	/* close files */
2861	SafeFclose (&fpParts);
2862	SafeFclose (&fpTstat);
2863	SafeFclose (&fpVstat);
2864	SafeFclose (&fpCon);
2865	SafeFclose (&fpTrees);
2866
2867	#if defined (PRINT_RATEMULTIPLIERS_CPP)
2868	SafeFclose (&rateMultfp);
2869	#endif
2870
2871	/* free pointer array to partitions */
2872	free (treeParts);
2873	treeParts = NULL;
2874	FreePartCtr (partCtrRoot);
2875	partCtrRoot = NULL;
2876	FreeTreeCtr (treeCtrRoot);
2877	treeCtrRoot = NULL;
2878	} /* next tree */
2879
2880	/* free memory and file pointers */
2881	if (s) free(s);
2882	FreeSumtParams();
2883
2884	/* reset numLocalTaxa and localOutGroup */
2885	ResetTaxonSet();
2886
2887	# if defined (MPI_ENABLED)
2888	}
2889	# endif
2890
2891	expecting = Expecting(COMMAND);
2892	inSumtCommand = NO;
2893	SafeFree ((void **)&tempStr);
2894
2895	return (NO_ERROR);
2896
2897	/* error exit */
2898	errorExit:
2899	/* free sumtParams */
2900	if (s) free(s);
2901	FreeSumtParams();
2902
2903	/* close files in case they are open*/
2904	SafeFclose (&fp);
2905	SafeFclose (&fpParts);
2906	SafeFclose (&fpTstat);
2907	SafeFclose (&fpVstat);
2908	SafeFclose (&fpCon);
2909	SafeFclose (&fpTrees);
2910
2911	#if defined (PRINT_RATEMULTIPLIERS_CPP)
2912	SafeFclose (&rateMultfp);
2913	#endif
2914
2915	/* free pointer array to partitions, part and tree counters */
2916	free (treeParts);
2917	FreePartCtr (partCtrRoot);
2918	FreeTreeCtr (treeCtrRoot);
2919	partCtrRoot = NULL;
2920	treeCtrRoot = NULL;
2921
2922	/* reset taxon set */
2923	ResetTaxonSet();
2924
2925	expecting = Expecting(COMMAND);
2926	inSumtCommand = NO;
2927	SafeFree ((void **)&tempStr);
2928
2929	return (ERROR);
2930	}
2931
2932
2933
2934
2935	void PrintSumtTableLine(int numRuns, int rowCount, Stat theStats, MrBFlt numPSRFSamples, MrBFlt maxPSRF, MrBFlt *sumPSRF)
2936	{
2937	int j,k;
2938
2939	MrBayesPrint ("%10.6lf %10.6lf %10.6lf %10.6lf %10.6lf", theStats->mean, theStats->var, theStats->lower, theStats->upper, theStats->median);
2940
2941	MrBayesPrintf (fpVstat, "\t%s", MbPrintNum(theStats->mean));
2942	MrBayesPrintf (fpVstat, "\t%s", MbPrintNum(theStats->var));
2943	MrBayesPrintf (fpVstat, "\t%s", MbPrintNum(theStats->lower));
2944	MrBayesPrintf (fpVstat, "\t%s", MbPrintNum(theStats->upper));
2945	MrBayesPrintf (fpVstat, "\t%s", MbPrintNum(theStats->median));
2946
2947	if (numRuns > 1)
2948	{
2949	for (j=k=0; j<numRuns; j++)
2950	if (rowCount[j] > 0)
2951	k++;
2952	if (theStats->PSRF < 0.0)
2953	{
2954	MrBayesPrint (" NA %3d", k);
2955	MrBayesPrintf (fpVstat, "\tNA\t%d", k);
2956	}
2957	else
2958	{
2959	if(theStats->PSRF > 10.0)
2960	{
2961	MrBayesPrint (" >10.0 %3d", k);
2962	MrBayesPrintf (fpVstat, "\tNA\t%d", k);
2963	(*maxPSRF) = 10.0;
2964	}
2965	else
2966	{
2967	MrBayesPrint (" %7.3lf %3d", theStats->PSRF, k);
2968	MrBayesPrintf (fpVstat, "\t%s\t%d", MbPrintNum(theStats->PSRF), k);
2969	(*sumPSRF) += theStats->PSRF;
2970	(*numPSRFSamples)++;
2971	if (theStats->PSRF > *maxPSRF)
2972	(*maxPSRF) = theStats->PSRF;
2973	}
2974	}
2975
2976	if (k != numRuns)
2977	MrBayesPrint (" *");
2978	}
2979
2980	MrBayesPrintf (fpVstat, "\n");
2981	MrBayesPrint ("\n");
2982	}
2983
2984
2985	int DoSumtParm (char parmName, char tkn)
2986
2987	{
2988
2989	int tempI;
2990	MrBFlt tempD;
2991	char tempStr[100];
2992
2993	if (defMatrix == NO)
2994	{
2995	MrBayesPrint ("%s A matrix must be specified before sumt can be used\n", spacer);
2996	return (ERROR);
2997	}
2998
2999	if (expecting == Expecting(PARAMETER))
3000	{
3001	expecting = Expecting(EQUALSIGN);
3002	}
3003	else
3004	{
3005	if (!strcmp(parmName, "Xxxxxxxxxx"))
3006	{
3007	expecting = Expecting(PARAMETER);
3008	expecting \|= Expecting(SEMICOLON);
3009	}
3010	/* set Filename (sumtParams.sumtFileName) ***************************************************/
3011	else if (!strcmp(parmName, "Filename"))
3012	{
3013	if (expecting == Expecting(EQUALSIGN))
3014	{
3015	expecting = Expecting(ALPHA);
3016	readWord = YES;
3017	}
3018	else if (expecting == Expecting(ALPHA))
3019	{
3020	if(strlen(tkn)>99)
3021	{
3022	MrBayesPrint ("%s Maximum allowed length of file name is 99 characters. The given name:\n", spacer);
3023	MrBayesPrint ("%s '%s'\n", spacer,tkn);
3024	MrBayesPrint ("%s has %d characters.\n", spacer,strlen(tkn));
3025	return (ERROR);
3026	}
3027	strcpy (sumtParams.sumtFileName, tkn);
3028	strcpy(sumtParams.sumtOutfile, tkn);
3029	MrBayesPrint ("%s Setting sumt filename and outputname to %s\n", spacer, sumtParams.sumtFileName);
3030	expecting = Expecting(PARAMETER) \| Expecting(SEMICOLON);
3031	}
3032	else
3033	return (ERROR);
3034	}
3035	/* set Relburnin (chainParams.relativeBurnin) ********************************************************/
3036	else if (!strcmp(parmName, "Relburnin"))
3037	{
3038	if (expecting == Expecting(EQUALSIGN))
3039	expecting = Expecting(ALPHA);
3040	else if (expecting == Expecting(ALPHA))
3041	{
3042	if (IsArgValid(tkn, tempStr) == NO_ERROR)
3043	{
3044	if (!strcmp(tempStr, "Yes"))
3045	chainParams.relativeBurnin = YES;
3046	else
3047	chainParams.relativeBurnin = NO;
3048	}
3049	else
3050	{
3051	MrBayesPrint ("%s Invalid argument for Relburnin\n", spacer);
3052	return (ERROR);
3053	}
3054	if (chainParams.relativeBurnin == YES)
3055	MrBayesPrint ("%s Using relative burnin (a fraction of samples discarded).\n", spacer);
3056	else
3057	MrBayesPrint ("%s Using absolute burnin (a fixed number of samples discarded).\n", spacer);
3058	expecting = Expecting(PARAMETER) \| Expecting(SEMICOLON);
3059	}
3060	else
3061	{
3062	return (ERROR);
3063	}
3064	}
3065	/* set Burnin (chainParams.chainBurnIn) ***********************************************************/
3066	else if (!strcmp(parmName, "Burnin"))
3067	{
3068	if (expecting == Expecting(EQUALSIGN))
3069	expecting = Expecting(NUMBER);
3070	else if (expecting == Expecting(NUMBER))
3071	{
3072	sscanf (tkn, "%d", &tempI);
3073	chainParams.chainBurnIn = tempI;
3074	MrBayesPrint ("%s Setting urn-in to %d\n", spacer, chainParams.chainBurnIn);
3075	expecting = Expecting(PARAMETER) \| Expecting(SEMICOLON);
3076	}
3077	else
3078	{
3079	return (ERROR);
3080	}
3081	}
3082	/* set Burninfrac (chainParams.burninFraction) ************************************************************/
3083	else if (!strcmp(parmName, "Burninfrac"))
3084	{
3085	if (expecting == Expecting(EQUALSIGN))
3086	expecting = Expecting(NUMBER);
3087	else if (expecting == Expecting(NUMBER))
3088	{
3089	sscanf (tkn, "%lf", &tempD);
3090	if (tempD < 0.01)
3091	{
3092	MrBayesPrint ("%s Burnin fraction too low (< 0.01)\n", spacer);
3093	return (ERROR);
3094	}
3095	if (tempD > 0.50)
3096	{
3097	MrBayesPrint ("%s Burnin fraction too high (> 0.50)\n", spacer);
3098	return (ERROR);
3099	}
3100	chainParams.burninFraction = tempD;
3101	MrBayesPrint ("%s Setting burnin fraction to %.2f\n", spacer, chainParams.burninFraction);
3102	expecting = Expecting(PARAMETER) \| Expecting(SEMICOLON);
3103	}
3104	else
3105	{
3106	return (ERROR);
3107	}
3108	}
3109	/* set Nruns (sumtParams.numRuns) *******************************************************/
3110	else if (!strcmp(parmName, "Nruns"))
3111	{
3112	if (expecting == Expecting(EQUALSIGN))
3113	expecting = Expecting(NUMBER);
3114	else if (expecting == Expecting(NUMBER))
3115	{
3116	sscanf (tkn, "%d", &tempI);
3117	if (tempI < 1)
3118	{
3119	MrBayesPrint ("%s Nruns must be at least 1\n", spacer);
3120	return (ERROR);
3121	}
3122	else
3123	{
3124	sumtParams.numRuns = tempI;
3125	MrBayesPrint ("%s Setting sumt nruns to %ld\n", spacer, sumtParams.numRuns);
3126	expecting = Expecting(PARAMETER) \| Expecting(SEMICOLON);
3127	}
3128	}
3129	else
3130	return (ERROR);
3131	}
3132	/* set Ntrees (sumtParams.numTrees) *******************************************************/
3133	else if (!strcmp(parmName, "Ntrees"))
3134	{
3135	if (expecting == Expecting(EQUALSIGN))
3136	expecting = Expecting(NUMBER);
3137	else if (expecting == Expecting(NUMBER))
3138	{
3139	sscanf (tkn, "%d", &tempI);
3140	if (tempI < 1)
3141	{
3142	MrBayesPrint ("%s Ntrees must be at least 1\n", spacer);
3143	return (ERROR);
3144	}
3145	else
3146	{
3147	sumtParams.numTrees = tempI;
3148	MrBayesPrint ("%s Setting sumt ntrees to %ld\n", spacer, sumtParams.numTrees);
3149	expecting = Expecting(PARAMETER) \| Expecting(SEMICOLON);
3150	}
3151	}
3152	else
3153	return (ERROR);
3154	}
3155	/* set Contype (sumtParams.sumtConType) *****************************************************/
3156	else if (!strcmp(parmName, "Contype"))
3157	{
3158	if (expecting == Expecting(EQUALSIGN))
3159	expecting = Expecting(ALPHA);
3160	else if (expecting == Expecting(ALPHA))
3161	{
3162	if (IsArgValid(tkn, tempStr) == NO_ERROR)
3163	{
3164	strcpy (sumtParams.sumtConType, tempStr);
3165	MrBayesPrint ("%s Setting sumt contype to %s\n", spacer, sumtParams.sumtConType);
3166	}
3167	expecting = Expecting(PARAMETER) \| Expecting(SEMICOLON);
3168	}
3169	else
3170	return (ERROR);
3171	}
3172	/* set Conformat (sumtParams.consensusFormat) *****************************************************/
3173	else if (!strcmp(parmName, "Conformat"))
3174	{
3175	if (expecting == Expecting(EQUALSIGN))
3176	expecting = Expecting(ALPHA);
3177	else if (expecting == Expecting(ALPHA))
3178	{
3179	if (IsArgValid(tkn, tempStr) == NO_ERROR)
3180	{
3181	if (!strcmp(tempStr,"Figtree"))
3182	{
3183	MrBayesPrint ("%s Setting sumt conformat to Figtree\n", spacer);
3184	sumtParams.consensusFormat = FIGTREE;
3185	}
3186	else
3187	{
3188	MrBayesPrint ("%s Setting sumt conformat to Simple\n", spacer);
3189	sumtParams.consensusFormat = SIMPLE;
3190	}
3191	}
3192	else
3193	{
3194	MrBayesPrint ("%s Invalid argument for calctreeprobs\n", spacer);
3195	return (ERROR);
3196	}
3197	expecting = Expecting(PARAMETER) \| Expecting(SEMICOLON);
3198	}
3199	else
3200	return (ERROR);
3201	}
3202	/* set Calctreeprobs (sumtParams.calcTreeprobs) *********************************************/
3203	else if (!strcmp(parmName, "Calctreeprobs"))
3204	{
3205	if (expecting == Expecting(EQUALSIGN))
3206	expecting = Expecting(ALPHA);
3207	else if (expecting == Expecting(ALPHA))
3208	{
3209	if (IsArgValid(tkn, tempStr) == NO_ERROR)
3210	{
3211	if (!strcmp(tempStr, "Yes"))
3212	sumtParams.calcTreeprobs = YES;
3213	else
3214	sumtParams.calcTreeprobs = NO;
3215	}
3216	else
3217	{
3218	MrBayesPrint ("%s Invalid argument for calctreeprobs\n", spacer);
3219	return (ERROR);
3220	}
3221	if (sumtParams.calcTreeprobs == YES)
3222	MrBayesPrint ("%s Setting calctreeprobs to yes\n", spacer);
3223	else
3224	MrBayesPrint ("%s Setting calctreeprobs to no\n", spacer);
3225	expecting = Expecting(PARAMETER) \| Expecting(SEMICOLON);
3226	}
3227	else
3228	return (ERROR);
3229	}
3230	/* set Showtreeprobs (sumtParams.showSumtTrees) *********************************************/
3231	else if (!strcmp(parmName, "Showtreeprobs"))
3232	{
3233	if (expecting == Expecting(EQUALSIGN))
3234	expecting = Expecting(ALPHA);
3235	else if (expecting == Expecting(ALPHA))
3236	{
3237	if (IsArgValid(tkn, tempStr) == NO_ERROR)
3238	{
3239	if (!strcmp(tempStr, "Yes"))
3240	sumtParams.showSumtTrees = YES;
3241	else
3242	sumtParams.showSumtTrees = NO;
3243	}
3244	else
3245	{
3246	MrBayesPrint ("%s Invalid argument for showtreeprobs\n", spacer);
3247	return (ERROR);
3248	}
3249	if (sumtParams.showSumtTrees == YES)
3250	MrBayesPrint ("%s Setting showtreeprobs to yes\n", spacer);
3251	else
3252	MrBayesPrint ("%s Setting showtreeprobs to no\n", spacer);
3253	expecting = Expecting(PARAMETER) \| Expecting(SEMICOLON);
3254	}
3255	else
3256	return (ERROR);
3257	}
3258	/* set Hpd (sumpParams.HPD) ********************************************************/
3259	else if (!strcmp(parmName, "Hpd"))
3260	{
3261	if (expecting == Expecting(EQUALSIGN))
3262	expecting = Expecting(ALPHA);
3263	else if (expecting == Expecting(ALPHA))
3264	{
3265	if (IsArgValid(tkn, tempStr) == NO_ERROR)
3266	{
3267	if (!strcmp(tempStr, "Yes"))
3268	sumtParams.HPD = YES;
3269	else
3270	sumtParams.HPD = NO;
3271	}
3272	else
3273	{
3274	MrBayesPrint ("%s Invalid argument for Hpd\n", spacer);
3275	return (ERROR);
3276	}
3277	if (sumtParams.HPD == YES)
3278	MrBayesPrint ("%s Reporting 95 %% region of Highest Posterior Density (HPD).\n", spacer);
3279	else
3280	MrBayesPrint ("%s Reporting median interval containing 95 %% of values.\n", spacer);
3281	expecting = Expecting(PARAMETER) \| Expecting(SEMICOLON);
3282	}
3283	else
3284	return (ERROR);
3285	}
3286	/* set Printbrlens (sumtParams.printBrlensToFile) *********************************************/
3287	else if (!strcmp(parmName, "Printbrlens"))
3288	{
3289	if (expecting == Expecting(EQUALSIGN))
3290	expecting = Expecting(ALPHA);
3291	else if (expecting == Expecting(ALPHA))
3292	{
3293	if (IsArgValid(tkn, tempStr) == NO_ERROR)
3294	{
3295	if (!strcmp(tempStr, "Yes"))
3296	sumtParams.printBrlensToFile = YES;
3297	else
3298	sumtParams.printBrlensToFile = NO;
3299	}
3300	else
3301	{
3302	MrBayesPrint ("%s Invalid argument for printbrlens\n", spacer);
3303	return (ERROR);
3304	}
3305	if (sumtParams.printBrlensToFile == YES)
3306	MrBayesPrint ("%s Setting printbrlens to yes\n", spacer);
3307	else
3308	MrBayesPrint ("%s Setting printbrlens to no\n", spacer);
3309	expecting = Expecting(PARAMETER) \| Expecting(SEMICOLON);
3310	}
3311	else
3312	return (ERROR);
3313	}
3314	/* set Brlensgeq (sumtParams.brlensFreqDisplay) *******************************************************/
3315	else if (!strcmp(parmName, "Brlensgeq"))
3316	{
3317	if (expecting == Expecting(EQUALSIGN))
3318	expecting = Expecting(NUMBER);
3319	else if (expecting == Expecting(NUMBER))
3320	{
3321	sscanf (tkn, "%lf", &tempD);
3322	sumtParams.brlensFreqDisplay = tempD;
3323	MrBayesPrint ("%s Printing branch lengths to file for partitions with probability >= %lf\n", spacer, sumtParams.brlensFreqDisplay);
3324	expecting = Expecting(PARAMETER) \| Expecting(SEMICOLON);
3325	}
3326	else
3327	return (ERROR);
3328	}
3329	/* set Ordertaxa (sumtParams.orderTaxa) *********************************************/
3330	else if (!strcmp(parmName, "Ordertaxa"))
3331	{
3332	if (expecting == Expecting(EQUALSIGN))
3333	expecting = Expecting(ALPHA);
3334	else if (expecting == Expecting(ALPHA))
3335	{
3336	if (IsArgValid(tkn, tempStr) == NO_ERROR)
3337	{
3338	if (!strcmp(tempStr, "Yes"))
3339	sumtParams.orderTaxa = YES;
3340	else
3341	sumtParams.orderTaxa = NO;
3342	}
3343	else
3344	{
3345	MrBayesPrint ("%s Invalid argument for ordertaxa\n", spacer);
3346	return (ERROR);
3347	}
3348	if (sumtParams.orderTaxa == YES)
3349	MrBayesPrint ("%s Setting ordertaxa to yes\n", spacer);
3350	else
3351	MrBayesPrint ("%s Setting ordertaxa to no\n", spacer);
3352	expecting = Expecting(PARAMETER) \| Expecting(SEMICOLON);
3353	}
3354	else
3355	return (ERROR);
3356	}
3357	/* set Outputname (sumtParams.sumtOutfile) *******************************************************/
3358	else if (!strcmp(parmName, "Outputname"))
3359	{
3360	if (expecting == Expecting(EQUALSIGN))
3361	{
3362	expecting = Expecting(ALPHA);
3363	readWord = YES;
3364	}
3365	else if (expecting == Expecting(ALPHA))
3366	{
3367	sscanf (tkn, "%s", tempStr);
3368	strcpy (sumtParams.sumtOutfile, tempStr);
3369	MrBayesPrint ("%s Setting sumt output file name to \"%s\"\n", spacer, sumtParams.sumtOutfile);
3370	expecting = Expecting(PARAMETER) \| Expecting(SEMICOLON);
3371	}
3372	else
3373	return (ERROR);
3374	}
3375	/* set Table (sumtParams.table) ********************************************************/
3376	else if (!strcmp(parmName, "Table"))
3377	{
3378	if (expecting == Expecting(EQUALSIGN))
3379	expecting = Expecting(ALPHA);
3380	else if (expecting == Expecting(ALPHA))
3381	{
3382	if (IsArgValid(tkn, tempStr) == NO_ERROR)
3383	{
3384	if (!strcmp(tempStr, "Yes"))
3385	sumtParams.table = YES;
3386	else
3387	sumtParams.table = NO;
3388	}
3389	else
3390	{
3391	MrBayesPrint ("%s Invalid argument for Table (valid arguments are 'yes' and 'no')\n", spacer);
3392	return (ERROR);
3393	}
3394	if (sumtParams.table == YES)
3395	MrBayesPrint ("%s Setting sumt to compute table of partition frequencies\n", spacer);
3396	else
3397	MrBayesPrint ("%s Setting sumt not to compute table of partition frequencies\n", spacer);
3398	expecting = Expecting(PARAMETER) \| Expecting(SEMICOLON);
3399	}
3400	else
3401	return (ERROR);
3402	}
3403	/* set Summary (sumtParams.summary) ********************************************************/
3404	else if (!strcmp(parmName, "Summary"))
3405	{
3406	if (expecting == Expecting(EQUALSIGN))
3407	expecting = Expecting(ALPHA);
3408	else if (expecting == Expecting(ALPHA))
3409	{
3410	if (IsArgValid(tkn, tempStr) == NO_ERROR)
3411	{
3412	if (!strcmp(tempStr, "Yes"))
3413	sumtParams.summary = YES;
3414	else
3415	sumtParams.summary = NO;
3416	}
3417	else
3418	{
3419	MrBayesPrint ("%s Invalid argument for 'Summary' (valid arguments are 'yes' and 'no')\n", spacer);
3420	return (ERROR);
3421	}
3422	if (sumtParams.summary == YES)
3423	MrBayesPrint ("%s Setting sumt to summary statistics\n", spacer);
3424	else
3425	MrBayesPrint ("%s Setting sumt not to compute summary statistics\n", spacer);
3426	expecting = Expecting(PARAMETER) \| Expecting(SEMICOLON);
3427	}
3428	else
3429	return (ERROR);
3430	}
3431	/* set Consensus (sumtParams.showConsensus) ********************************************************/
3432	else if (!strcmp(parmName, "Consensus"))
3433	{
3434	if (expecting == Expecting(EQUALSIGN))
3435	expecting = Expecting(ALPHA);
3436	else if (expecting == Expecting(ALPHA))
3437	{
3438	if (IsArgValid(tkn, tempStr) == NO_ERROR)
3439	{
3440	if (!strcmp(tempStr, "Yes"))
3441	sumtParams.showConsensus = YES;
3442	else
3443	sumtParams.showConsensus = NO;
3444	}
3445	else
3446	{
3447	MrBayesPrint ("%s Invalid argument for Consensus (valid arguments are 'yes' and 'no')\n", spacer);
3448	return (ERROR);
3449	}
3450	if (sumtParams.showConsensus == YES)
3451	MrBayesPrint ("%s Setting sumt to show consensus trees\n", spacer);
3452	else
3453	MrBayesPrint ("%s Setting sumt not to show consensus trees\n", spacer);
3454	expecting = Expecting(PARAMETER) \| Expecting(SEMICOLON);
3455	}
3456	else
3457	return (ERROR);
3458	}
3459	/* set Minpartfreq (sumtParams.minPartFreq) *******************************************************/
3460	else if (!strcmp(parmName, "Minpartfreq"))
3461	{
3462	if (expecting == Expecting(EQUALSIGN))
3463	expecting = Expecting(NUMBER);
3464	else if (expecting == Expecting(NUMBER))
3465	{
3466	sscanf (tkn, "%lf", &tempD);
3467	sumtParams.minPartFreq = tempD;
3468	MrBayesPrint ("%s Including partitions with probability greater than or equal to %lf in summary statistics\n", spacer, sumtParams.minPartFreq);
3469	expecting = Expecting(PARAMETER) \| Expecting(SEMICOLON);
3470	}
3471	else
3472	return (ERROR);
3473	}
3474	else
3475	return (ERROR);
3476	}
3477
3478	return (NO_ERROR);
3479
3480	}
3481
3482
3483
3484
3485
3486	int DoSumtTree (void)
3487
3488	{
3489
3490	int i, j, z, printEvery, nAstPerPrint, burnin;
3491	MrBFlt x, y;
3492	PolyTree *t;
3493	PolyNode *p;
3494
3495	#if defined (PRINT_RATEMULTIPLIERS_CPP)
3496
3497	/* get depths if relevant */
3498	if (sumtParams.tree->isClock)
3499	GetPolyDepths (sumtParams.tree);
3500
3501	if(rateMultfp!=NULL && sumtParams.tree->root!=NULL)
3502	DELETE_ME_dump_depth(sumtParams.tree->root);
3503	//fprintf(rateMultfp,"%s\n",tkn);
3504	#endif
3505
3506	/* increment number of trees read in */
3507	sumtParams.numFileTrees[sumtParams.runId]++;
3508	sumtParams.numTreesEncountered++;
3509
3510	/* update status bar */
3511	if (sumtParams.numTreesInLastBlock * sumtParams.numRuns < 80)
3512	{
3513	printEvery = 1;
3514	nAstPerPrint = 80 / (sumtParams.numTreesInLastBlock * sumtParams.numRuns);
3515	if (sumtParams.numTreesEncountered % printEvery == 0)
3516	{
3517	for (i=0; i<nAstPerPrint; i++)
3518	{
3519	MrBayesPrint ("*");
3520	numAsterices++;
3521	}
3522	}
3523	}
3524	else
3525	{
3526	x = (MrBFlt)(sumtParams.numTreesInLastBlock * sumtParams.numRuns) / (MrBFlt) (80);
3527	y = (MrBFlt)(sumtParams.numFileTrees[sumtParams.runId] + sumtParams.numTreesInLastBlock * sumtParams.runId) / x;
3528	z = (int)y;
3529	if (numAsterices < z)
3530	{
3531	MrBayesPrint ("*");
3532	numAsterices++;
3533	}
3534	}
3535
3536	/* get burnin */
3537	if (inComparetreeCommand == YES)
3538	burnin = comptreeParams.burnin;
3539	else
3540	burnin = sumtParams.burnin;
3541
3542	if (sumtParams.numFileTrees[sumtParams.runId] > burnin)
3543	{
3544	/* increment the number of trees sampled */
3545	sumtParams.numFileTreesSampled[sumtParams.runId]++;
3546	sumtParams.numTreesSampled++;
3547
3548	/* get the tree we just read in */
3549	t = sumtParams.tree;
3550
3551	/* move calculation root for nonrooted trees if necessary */
3552	MovePolyCalculationRoot (t, localOutGroup);
3553
3554	/* check taxon set and outgroup */
3555	if (sumtParams.runId == 0 && sumtParams.numFileTreesSampled[0] == 1)
3556	{
3557	if (isTranslateDef == YES && isTranslateDiff == YES)
3558	{
3559	/* we are using a translate block with different taxa set */
3560	if (t->nNodes - t->nIntNodes != numTranslates)
3561	{
3562	MrBayesPrint ("%s ERROR: Expected %d taxa; found %d taxa\n", spacer, numTranslates, sumtParams.numTaxa);
3563	return (ERROR);
3564	}
3565	for (i=0; i<numTaxa; i++)
3566	sumtParams.absentTaxa[i] = NO;
3567	localOutGroup = 0; /* no previous outgroup assignment is valid */
3568	}
3569	else
3570	{
3571	/* we are using the current taxa set */
3572	for (i=0; i<numTaxa; i++)
3573	sumtParams.absentTaxa[i] = YES;
3574	for (i=0; i<t->nNodes; i++)
3575	{
3576	p = t->allDownPass[i];
3577	if (p->left == NULL)
3578	sumtParams.absentTaxa[p->index] = NO;
3579	}
3580	localOutGroup = 0;
3581	for (i=j=0; i<numTaxa; i++)
3582	{
3583	if (sumtParams.absentTaxa[i] == NO && taxaInfo[i].isDeleted == NO)
3584	{
3585	if (i == outGroupNum)
3586	localOutGroup = j;
3587	j++;
3588	}
3589	}
3590	}
3591
3592	/* now we can safely prune the tree based on taxaInfo[].isDeleted */
3593	/* the following block was conditioned with if(isTranslateDef == NO \|\| isTranslateDiff == NO)
3594	The reason was not clearly stated but it prevents exclusion of taxa to work in case when the condition does not hold.
3595	My guess is that before PrunePolyTree() relied on indeses of tips be set as in original matrix.
3596	Now it is not needed after PrunePolyTree and ResetTipIndices ware modified to use labels istead of indexes to recognize tips.*/
3597	{
3598	PrunePolyTree(t);
3599
3600	/* reset tip and int node indices in case some taxa deleted */
3601	ResetTipIndices (t);
3602	ResetIntNodeIndices(t);
3603	}
3604
3605	/* set basic parameters */
3606	sumtParams.numTaxa = t->nNodes - t->nIntNodes;
3607	numLocalTaxa = sumtParams.numTaxa;
3608	sumtParams.SafeLongsNeeded = ((numLocalTaxa-1) / nBitsInALong) + 1;
3609	if (t->isRooted == YES)
3610	sumtParams.orderLen = numLocalTaxa - 2;
3611	else
3612	sumtParams.orderLen = numLocalTaxa - 3;
3613	}
3614	else
3615	{
3616	/* the following block was conditioned with if(isTranslateDef == NO \|\| isTranslateDiff == NO)
3617	The reason was not clearly stated but it prevents exclusion of taxa to work in case when the condition does not hold.
3618	My guess is that before PrunePolyTree() relied on indeses of tips be set as in original matrix.
3619	Now it is not needed after PrunePolyTree and ResetTipIndices ware modified to use labels istead of indexes to recognize tips.*/
3620	{
3621	for (i=0; i<t->nNodes; i++)
3622	{
3623	p = t->allDownPass[i];
3624	if (p->left == NULL && taxaInfo[p->index].isDeleted == NO && sumtParams.absentTaxa[p->index] == YES)
3625	{
3626	MrBayesPrint ("%s Taxon %d should not be in sampled tree\n", spacer, p->index + 1);
3627	return (ERROR);
3628	}
3629	}
3630
3631	/* now we can safely prune the tree based on taxaInfo[].isDeleted */
3632	PrunePolyTree (t);
3633
3634	/* reset tip and int node indices in case some taxa deleted */
3635	ResetTipIndices (t);
3636	ResetIntNodeIndices(t);
3637	}
3638
3639	/* check that all taxa are included */
3640	if (t->nNodes - t->nIntNodes != sumtParams.numTaxa)
3641	{
3642	MrBayesPrint ("%s Expecting %d taxa but tree '%s' in file '%s' has %d taxa\n",
3643	spacer, sumtParams.numTaxa, t->name, sumtParams.curFileName, t->nNodes-t->nIntNodes);
3644	return ERROR;
3645	}
3646	}
3647
3648	if (sumtParams.runId == 0 && sumtParams.numFileTreesSampled[0] == 1)
3649	{
3650	/* harvest labels (can only be done safely after pruning) */
3651	for (i=0; i<sumtParams.numTaxa; i++)
3652	{
3653	for (j=0; j<t->nNodes; j++)
3654	{
3655	p = t->allDownPass[j];
3656	if (p->index == i){
3657	if ( strlen(p->label)>99 )
3658	{
3659	MrBayesPrint ("%s Taxon name %s is too long. Maximun 99 characters is allowed.\n", spacer, p->label);
3660	return (ERROR);
3661	}
3662	AddString(&sumtParams.taxaNames, i, p->label);
3663	}
3664	}
3665	}
3666	}
3667
3668	/* check that tree agrees with template */
3669	if (sumtParams.numTreesSampled == 1)
3670	{
3671	sumtParams.brlensDef = t->brlensDef;
3672	sumtParams.isRooted = t->isRooted;
3673	sumtParams.isClock = t->isClock;
3674	sumtParams.isCalibrated = t->isCalibrated;
3675	sumtParams.isRelaxed = t->isRelaxed;
3676	sumtParams.nBSets = 0;
3677	sumtParams.nESets = 0;
3678	for (i=0; i<t->nBSets; i++)
3679	AddString(&sumtParams.bSetName,sumtParams.nBSets++,t->bSetName[i]);
3680	for (i=0; i<t->nESets; i++)
3681	AddString(&sumtParams.eSetName,sumtParams.nESets++,t->eSetName[i]);
3682	if (t->popSizeSet == YES)
3683	{
3684	sumtParams.popSizeSet = YES;
3685	sumtParams.popSizeSetName = (char *) SafeCalloc (strlen(t->popSizeSetName)+1, sizeof(char));
3686	strcpy(sumtParams.popSizeSetName, t->popSizeSetName);
3687	}
3688	else
3689	sumtParams.popSizeSet = NO;
3690	}
3691	else /* if (sumtParams.numTreesSampled > 1) */
3692	{
3693	if (sumtParams.brlensDef != t->brlensDef)
3694	{
3695	MrBayesPrint ("%s Trees with and without branch lengths mixed\n", spacer);
3696	return ERROR;
3697	}
3698	if (sumtParams.isRooted != t->isRooted)
3699	{
3700	if (sumtParams.isRooted == YES)
3701	MrBayesPrint ("%s Expected rooted tree but tree '%s' in file '%s' is not rooted\n",
3702	spacer, t->name, sumtParams.curFileName);
3703	else if (sumtParams.isRooted == NO)
3704	MrBayesPrint ("%s Expected unrooted tree but tree '%s' in file '%s' is rooted\n",
3705	spacer, t->name, sumtParams.curFileName);
3706	return ERROR;
3707	}
3708	if (sumtParams.isClock != t->isClock)
3709	{
3710	if (sumtParams.isClock == YES)
3711	MrBayesPrint ("%s Expected clock tree but tree '%s' in file '%s' is not clock\n",
3712	spacer, t->name, sumtParams.curFileName);
3713	else if (sumtParams.isClock == NO)
3714	MrBayesPrint ("%s Expected nonclock tree but tree '%s' in file '%s' is clock\n",
3715	spacer, t->name, sumtParams.curFileName);
3716	return ERROR;
3717	}
3718	if (sumtParams.isCalibrated != t->isCalibrated)
3719	{
3720	if (sumtParams.isCalibrated == YES)
3721	MrBayesPrint ("%s Expected calibrated tree but tree '%s' in file '%s' is not calibrated\n",
3722	spacer, t->name, sumtParams.curFileName);
3723	else if (sumtParams.isCalibrated == NO)
3724	MrBayesPrint ("%s Expected noncalibrated tree but tree '%s' in file '%s' is calibrated\n",
3725	spacer, t->name, sumtParams.curFileName);
3726	return ERROR;
3727	}
3728	if (inComparetreeCommand == NO && sumtParams.isRelaxed != t->isRelaxed)
3729	{
3730	if (sumtParams.isRelaxed == YES)
3731	MrBayesPrint ("%s Expected relaxed clock tree but tree '%s' in file '%s' is not relaxed\n",
3732	spacer, t->name, sumtParams.curFileName);
3733	else if (sumtParams.isRelaxed == NO)
3734	MrBayesPrint ("%s Expected unrooted tree but tree '%s' in file '%s' is rooted\n",
3735	spacer, t->name, sumtParams.curFileName);
3736	return ERROR;
3737	}
3738	if (inComparetreeCommand == NO && (sumtParams.nESets != t->nESets \|\| sumtParams.nBSets != t->nBSets) )
3739	{
3740	MrBayesPrint ("%s Tree '%s' in file '%s' does not have the expected relaxed clock parameters\n",
3741	spacer, t->name, sumtParams.curFileName);
3742	return ERROR;
3743	}
3744	}
3745
3746	/* set partitions for tree */
3747	ResetPolyTreePartitions(t);
3748
3749	/* get depths if relevant */
3750	if (t->isClock)
3751	GetPolyDepths (t);
3752
3753	/* get ages if relevant */
3754	if (t->isCalibrated)
3755	GetPolyAges (t);
3756
3757	/* add partitions to counters */
3758	for (i=0; i<t->nNodes; i++)
3759	{
3760	p = t->allDownPass[i];
3761	partCtrRoot = AddSumtPartition (partCtrRoot, t, p, sumtParams.runId);
3762	}
3763
3764	/* add the tree to relevant tree list */
3765	if (inSumtCommand == YES)
3766	{
3767	if (t->isRooted == YES)
3768	StoreRPolyTopology (t, sumtParams.order);
3769	else /* if (sumtParams.isRooted == NO) */
3770	StoreUPolyTopology (t, sumtParams.order);
3771	treeCtrRoot = AddSumtTree (treeCtrRoot, sumtParams.order);
3772	}
3773	else
3774	{
3775	i = sumtParams.numFileTreesSampled[sumtParams.runId] - 1;
3776	if (StoreSumtTree (packedTreeList[sumtParams.runId], i, t) == ERROR)
3777	return (ERROR);
3778	}
3779
3780	/* Display the tree nodes. */
3781	# if 0
3782	ShowPolyNodes(t);
3783	# endif
3784	}
3785
3786	return (NO_ERROR);
3787	}
3788
3789
3790
3791
3792
3793	int ExamineSumtFile (char fileName, SumtFileInfo sumtFileInfo, char treeName, int brlensDef)
3794	{
3795	int i, foundBegin, lineTerm, inTreeBlock, blockErrors, inSumtComment, lineNum, numTreesInBlock,
3796	tokenType;
3797	char sumtToken[100], s, sumtTokenP;
3798	FILE *fp;
3799
3800	/* open binary file */
3801	if ((fp = OpenBinaryFileR(fileName)) == NULL)
3802	return ERROR;
3803
3804	/* find out what type of line termination is used for file 1 */
3805	lineTerm = LineTermType (fp);
3806	if (lineTerm != LINETERM_MAC && lineTerm != LINETERM_DOS && lineTerm != LINETERM_UNIX)
3807	{
3808	MrBayesPrint ("%s Unknown line termination for file \"%s\"\n", spacer, fileName);
3809	return ERROR;
3810	}
3811
3812	/* find length of longest line in either file */
3813	sumtFileInfo->longestLineLength = LongestLine (fp);
3814	sumtFileInfo->longestLineLength += 10;
3815
3816	/* allocate a string long enough to hold a line */
3817	s = (char )SafeMalloc((size_t) (sumtFileInfo->longestLineLength sizeof(char)));
3818	if (!s)
3819	{
3820	MrBayesPrint ("%s Problem allocating string for examining file \"%s\"\n", spacer, fileName);
3821	return (ERROR);
3822	}
3823
3824	/* close binary file */
3825	SafeFclose (&fp);
3826
3827	foundBegin = inTreeBlock = blockErrors = inSumtComment = NO;
3828	lineNum = numTreesInBlock = 0;
3829	sumtFileInfo->numTreeBlocks = 0;
3830	sumtFileInfo->lastTreeBlockBegin = 0;
3831	sumtFileInfo->lastTreeBlockEnd = 0;
3832	sumtFileInfo->numTreesInLastBlock = 0;
3833
3834	/* open text file */
3835	if ((fp = OpenTextFileR(fileName))==NULL)
3836	{
3837	MrBayesPrint ("%s Could not read file \"%s\" in text mode \n", spacer, fileName);
3838	return (ERROR);
3839	}
3840
3841	/* read file */
3842	while (fgets (s, sumtFileInfo->longestLineLength-2, fp) != NULL)
3843	{
3844	sumtTokenP = &s[0];
3845	do
3846	{
3847	if(GetToken (sumtToken, &tokenType, &sumtTokenP))
3848	goto errorExit;
3849	if (IsSame("[", sumtToken) == SAME)
3850	inSumtComment = YES;
3851	if (IsSame("]", sumtToken) == SAME)
3852	inSumtComment = NO;
3853
3854	if (inSumtComment == YES)
3855	{
3856	if (IsSame ("Param", sumtToken) == SAME)
3857	{
3858	/* extract the tree name */
3859	if(GetToken (sumtToken, &tokenType, &sumtTokenP)) /* get the colon */
3860	goto errorExit;
3861	if(GetToken (sumtToken, &tokenType, &sumtTokenP)) /* get the tree name */
3862	goto errorExit;
3863	strcpy (treeName, sumtToken);
3864	if(GetToken (sumtToken, &tokenType, &sumtTokenP))
3865	goto errorExit;
3866	while (IsSame("]", sumtToken) != SAME)
3867	{
3868	strcat (treeName, sumtToken);
3869	if(GetToken (sumtToken, &tokenType, &sumtTokenP))
3870	goto errorExit;
3871	}
3872	inSumtComment = NO;
3873	}
3874	}
3875	else /* if (inSumtComment == NO) */
3876	{
3877	if (foundBegin == YES)
3878	{
3879	if (IsSame("Trees", sumtToken) == SAME)
3880	{
3881	numTreesInBlock = 0;
3882	inTreeBlock = YES;
3883	foundBegin = NO;
3884	sumtFileInfo->lastTreeBlockBegin = lineNum;
3885	}
3886	}
3887	else
3888	{
3889	if (IsSame("Begin", sumtToken) == SAME)
3890	{
3891	if (foundBegin == YES)
3892	{
3893	MrBayesPrint ("%s Found inappropriate \"Begin\" statement in file\n", spacer);
3894	blockErrors = YES;
3895	}
3896	foundBegin = YES;
3897	}
3898	else if (IsSame("End", sumtToken) == SAME)
3899	{
3900	if (inTreeBlock == YES)
3901	{
3902	sumtFileInfo->numTreeBlocks++;
3903	inTreeBlock = NO;
3904	sumtFileInfo->lastTreeBlockEnd = lineNum;
3905	}
3906	else
3907	{
3908	MrBayesPrint ("%s Found inappropriate \"End\" statement in file\n", spacer);
3909	blockErrors = YES;
3910	}
3911	sumtFileInfo->numTreesInLastBlock = numTreesInBlock;
3912	}
3913	else if (IsSame("Tree", sumtToken) == SAME)
3914	{
3915	if (inTreeBlock == YES)
3916	{
3917	numTreesInBlock++;
3918	if (numTreesInBlock == 1)
3919	{
3920	*brlensDef = NO;
3921	for (i=0; s[i]!='\0'; i++)
3922	{
3923	if (s[i] == ':')
3924	{
3925	*brlensDef = YES;
3926	break;
3927	}
3928	}
3929	}
3930	}
3931	else
3932	{
3933	MrBayesPrint ("%s Found a \"Tree\" statement that is not in a tree block\n", spacer);
3934	blockErrors = YES;
3935	}
3936	}
3937	}
3938	}
3939
3940	} while (*sumtToken);
3941	lineNum++;
3942	}
3943
3944	/* Now, check some aspects of the tree file, such as the number of tree blocks and whether they are properly terminated. */
3945	if (inTreeBlock == YES)
3946	{
3947	MrBayesPrint ("%s Unterminated tree block in file %s. You probably need to\n", spacer, fileName);
3948	MrBayesPrint ("%s add a new line to the end of the file with \"End;\" on it.\n", spacer);
3949	goto errorExit;
3950	}
3951	if (inSumtComment == YES)
3952	{
3953	MrBayesPrint ("%s Unterminated comment in file %s\n", spacer, fileName);
3954	goto errorExit;
3955	}
3956	if (blockErrors == YES)
3957	{
3958	MrBayesPrint ("%s Found formatting errors in file %s\n", spacer, fileName);
3959	goto errorExit;
3960	}
3961	if (sumtFileInfo->lastTreeBlockEnd < sumtFileInfo->lastTreeBlockBegin)
3962	{
3963	MrBayesPrint ("%s Problem reading tree file %s\n", spacer, fileName);
3964	goto errorExit;
3965	}
3966	if (sumtFileInfo->numTreesInLastBlock <= 0)
3967	{
3968	MrBayesPrint ("%s No trees were found in last tree block of file %s\n", spacer, fileName);
3969	goto errorExit;
3970	}
3971	free (s);
3972	return (NO_ERROR);
3973
3974	errorExit:
3975	free (s);
3976	return (ERROR);
3977	}
3978
3979
3980
3981
3982
3983	/* FreePartCtr: Recursively free partition counter nodes */
3984	void FreePartCtr (PartCtr *r)
3985
3986	{
3987	int i, j;
3988
3989	if (r==NULL)
3990	return;
3991
3992	FreePartCtr (r->left);
3993	FreePartCtr (r->right);
3994
3995	/* free relaxed clock parameters: eRate, nEvents, bRate */
3996	if (sumtParams.nESets > 0)
3997	{
3998	for (i=0; i<sumtParams.nESets; i++)
3999	{
4000	for (j=0; j<sumtParams.numRuns; j++)
4001	free (r->nEvents[i][j]);
4002	free (r->nEvents[i]);
4003	}
4004	free (r->nEvents);
4005	}
4006	if (sumtParams.nBSets > 0)
4007	{
4008	for (i=0; i<sumtParams.nBSets; i++)
4009	{
4010	for (j=0; j<sumtParams.numRuns; j++)
4011	{
4012	free (r->bLen [i][j]);
4013	free (r->bRate[i][j]);
4014	}
4015	free (r->bLen [i]);
4016	free (r->bRate[i]);
4017	}
4018	free (r->bLen );
4019	free (r->bRate);
4020	}
4021
4022	/* free basic parameters */
4023	for (i=0; i<sumtParams.numRuns; i++)
4024	free (r->length[i]);
4025
4026	free (r->length);
4027	free (r->count);
4028	free (r->partition);
4029	free (r);
4030	numUniqueSplitsFound--;
4031	r = NULL;
4032	}
4033
4034
4035
4036
4037
4038	/* FreeSumtParams: Free parameters allocated in sumtParams struct */
4039	void FreeSumtParams(void)
4040	{
4041	int i;
4042
4043	if (memAllocs[ALLOC_SUMTPARAMS] == YES)
4044	{
4045	for (i=0; i<sumtParams.numTaxa; i++)
4046	free(sumtParams.taxaNames[i]);
4047	free (sumtParams.taxaNames);
4048	sumtParams.taxaNames = NULL;
4049	if (sumtParams.numFileTrees) free (sumtParams.numFileTrees);
4050	sumtParams.numFileTrees = NULL;
4051	FreePolyTree (sumtParams.tree);
4052	sumtParams.tree = NULL;
4053	if (sumtParams.nBSets > 0)
4054	{
4055	for (i=0; i<sumtParams.nBSets; i++)
4056	free(sumtParams.bSetName[i]);
4057	free (sumtParams.bSetName);
4058	sumtParams.bSetName = NULL;
4059	sumtParams.nBSets = 0;
4060	}
4061	if (sumtParams.nESets > 0)
4062	{
4063	for (i=0; i<sumtParams.nESets; i++)
4064	free(sumtParams.eSetName[i]);
4065	free (sumtParams.eSetName);
4066	sumtParams.eSetName = NULL;
4067	sumtParams.nESets = 0;
4068	}
4069	if (sumtParams.popSizeSet == YES)
4070	{
4071	free (sumtParams.popSizeSetName);
4072	sumtParams.popSizeSetName = NULL;
4073	sumtParams.popSizeSet = NO;
4074	}
4075	memAllocs[ALLOC_SUMTPARAMS] = NO;
4076	}
4077	}
4078
4079
4080
4081
4082
4083	/* FreeTreeCtr: Recursively free tree counter nodes */
4084	void FreeTreeCtr (TreeCtr *r)
4085
4086	{
4087
4088	if (r==NULL)
4089	return;
4090
4091	FreeTreeCtr (r->left);
4092	FreeTreeCtr (r->right);
4093
4094	free (r->order);
4095	free (r);
4096	numUniqueTreesFound--;
4097	r = NULL;
4098	}
4099
4100
4101
4102
4103
4104	/* Label: Calculate length of label and fill in char label if not NULL /
4105	int Label (PolyNode p, int addIndex, char label, int maxLength)
4106	{
4107	int i, j0, j1, k, n, length, nameLength, index;
4108
4109	if (p == NULL)
4110	return 0;
4111
4112	/* first calculate length */
4113	if (inSumtCommand == YES && isTranslateDiff == NO)
4114	{
4115	for (index=i=0; index<numTaxa; index++)
4116	{
4117	if (sumtParams.absentTaxa[index] == YES \|\| taxaInfo[index].isDeleted == YES)
4118	continue;
4119	if (p->index == i)
4120	break;
4121	else
4122	i++;
4123	}
4124	}
4125	else
4126	index = p->index;
4127
4128	if (addIndex != NO)
4129	length = (int)(strlen(p->label)) + 4 + (int)(log10(index+1));
4130	else
4131	length = (int)(strlen(p->label));
4132	length = (length > maxLength ? maxLength : length);
4133
4134	/* fill in label if label != NULL */
4135	if (label != NULL)
4136	{
4137	if (addIndex != NO)
4138	nameLength = length - 4 - (int)(log10(index+1));
4139	else
4140	nameLength = length;
4141
4142	for (i=0; i<nameLength-1; i++)
4143	label[i] = p->label[i];
4144	if ((int)strlen(p->label) > nameLength)
4145	label[i] = '~';
4146	else
4147	label[i] = p->label[i];
4148
4149	if (addIndex != NO)
4150	{
4151	label[++i] = ' ';
4152	label[++i] = '(';
4153	n = index + 1;
4154	k = (int)(log10(n)) + 1;
4155	while (n != 0)
4156	{
4157	j0 = (int)(log10(n));
4158	j1 = (int)(pow(10,j0));
4159	label[++i] = '0' + n/j1;
4160	n = n % j1;
4161	k--;
4162	}
4163	while (k!=0)
4164	{
4165	label[++i] = '0';
4166	k--;
4167	}
4168	label[++i] = ')';
4169	}
4170	label[++i] = '\0';
4171	}
4172
4173	return length;
4174	}
4175
4176
4177
4178
4179
4180	int OpenComptFiles (void)
4181
4182	{
4183
4184	int len, previousFiles, oldNoWarn, oldAutoOverwrite;
4185	char pFilename[100], dFilename[100];
4186	FILE *fpTemp;
4187
4188	oldNoWarn = noWarn;
4189	oldAutoOverwrite = autoOverwrite;
4190
4191	/* set file names */
4192	strcpy (pFilename, comptreeParams.comptOutfile);
4193	strcpy (dFilename, comptreeParams.comptOutfile);
4194	strcat (pFilename, ".pairs");
4195	strcat (dFilename, ".dists");
4196
4197	/* one overwrite check for both files */
4198	previousFiles = NO;
4199	if (noWarn == NO)
4200	{
4201	if ((fpTemp = OpenTextFileR(pFilename)) != NULL)
4202	{
4203	previousFiles = YES;
4204	fclose(fpTemp);
4205	}
4206	if ((fpTemp = OpenTextFileR(dFilename)) != NULL)
4207	{
4208	previousFiles = YES;
4209	fclose(fpTemp);
4210	}
4211	if (previousFiles == YES)
4212	{
4213	MrBayesPrint("%s There are previous compare results saved using the same filenames.\n", spacer);
4214	if (WantTo("Do you want to overwrite these results") == YES)
4215	{
4216	MrBayesPrint("\n");
4217	noWarn = YES;
4218	autoOverwrite = YES;
4219	}
4220	else
4221	{
4222	MrBayesPrint("\n");
4223	MrBayesPrint("%s Please specify a different output file name before running the comparetree command.\n", spacer);
4224	MrBayesPrint("%s You can do that using 'comparetree outputfile=<name>'. You can also move or\n", spacer);
4225	MrBayesPrint("%s rename the old result files.\n", spacer);
4226	return ERROR;
4227	}
4228	}
4229	}
4230
4231	if ((fpParts = OpenNewMBPrintFile (pFilename)) == NULL)
4232	{
4233	noWarn = oldNoWarn;
4234	autoOverwrite = oldAutoOverwrite;
4235	return ERROR;
4236	}
4237	if ((fpDists = OpenNewMBPrintFile (dFilename)) == NULL)
4238	{
4239	noWarn = oldNoWarn;
4240	autoOverwrite = oldAutoOverwrite;
4241	return ERROR;
4242	}
4243
4244	/* Reset file flags */
4245	noWarn = oldNoWarn;
4246	autoOverwrite = oldAutoOverwrite;
4247
4248	/* print unique identifiers to each file */
4249	len = (int) strlen (stamp);
4250	if (len > 1)
4251	{
4252	fprintf (fpParts, "[ID: %s]\n", stamp);
4253	fprintf (fpDists, "[ID: %s]\n", stamp);
4254	}
4255
4256	return (NO_ERROR);
4257	}
4258
4259
4260
4261
4262
4263	int OpenSumtFiles (int treeNo)
4264
4265	{
4266
4267	int i, len, oldNoWarn, oldAutoOverwrite, previousFiles;
4268	char pFilename[100], sFilename[100], vFilename[100], cFilename[100], tFilename[100];
4269	FILE *fpTemp;
4270
4271	oldNoWarn = noWarn;
4272	oldAutoOverwrite = autoOverwrite;
4273
4274	/* one overwrite check for all files */
4275	if (noWarn == NO && treeNo == 0)
4276	{
4277	previousFiles = NO;
4278	for (i=0; i<sumtParams.numTrees; i++)
4279	{
4280	if (sumtParams.numTrees > 1)
4281	{
4282	sprintf (pFilename, "%s.tree%d.parts", sumtParams.sumtOutfile, i+1);
4283	sprintf (sFilename, "%s.tree%d.tstat", sumtParams.sumtOutfile, i+1);
4284	sprintf (vFilename, "%s.tree%d.vstat", sumtParams.sumtOutfile, i+1);
4285	sprintf (cFilename, "%s.tree%d.con.tre", sumtParams.sumtOutfile, i+1);
4286	sprintf (tFilename, "%s.tree%d.trprobs", sumtParams.sumtOutfile, i+1);
4287	}
4288	else
4289	{
4290	sprintf (pFilename, "%s.parts", sumtParams.sumtOutfile);
4291	sprintf (sFilename, "%s.tstat", sumtParams.sumtOutfile);
4292	sprintf (vFilename, "%s.vstat", sumtParams.sumtOutfile);
4293	sprintf (cFilename, "%s.con.tre", sumtParams.sumtOutfile);
4294	sprintf (tFilename, "%s.trprobs", sumtParams.sumtOutfile);
4295	}
4296	if ((fpTemp = TestOpenTextFileR(pFilename)) != NULL)
4297	{
4298	previousFiles = YES;
4299	fclose(fpTemp);
4300	}
4301	if ((fpTemp = TestOpenTextFileR(sFilename)) != NULL)
4302	{
4303	previousFiles = YES;
4304	fclose(fpTemp);
4305	}
4306	if ((fpTemp = TestOpenTextFileR(vFilename)) != NULL)
4307	{
4308	previousFiles = YES;
4309	fclose(fpTemp);
4310	}
4311	if ((fpTemp = TestOpenTextFileR(cFilename)) != NULL)
4312	{
4313	previousFiles = YES;
4314	fclose(fpTemp);
4315	}
4316	if ((fpTemp = TestOpenTextFileR(tFilename)) != NULL)
4317	{
4318	previousFiles = YES;
4319	fclose(fpTemp);
4320	}
4321	if (previousFiles == YES)
4322	{
4323	MrBayesPrint("\n");
4324	MrBayesPrint("%s There are previous tree sample summaries saved using the same filenames.\n", spacer);
4325	if (WantTo("Do you want to overwrite these results") == YES)
4326	{
4327	MrBayesPrint("\n");
4328	noWarn = YES;
4329	autoOverwrite = YES;
4330	}
4331	else
4332	{
4333	MrBayesPrint("\n");
4334	MrBayesPrint("%s Please specify a different output file name before running the sumt command.\n", spacer);
4335	MrBayesPrint("%s You can do that using 'sumt outputfile=<name>'. You can also move or\n", spacer);
4336	MrBayesPrint("%s rename the old result files.\n", spacer);
4337	return ERROR;
4338	}
4339	}
4340	}
4341	}
4342
4343	/* set file names */
4344	if (sumtParams.numTrees > 1)
4345	{
4346	sprintf (pFilename, "%s.tree%d.parts", sumtParams.sumtOutfile, treeNo+1);
4347	sprintf (sFilename, "%s.tree%d.tstat", sumtParams.sumtOutfile, treeNo+1);
4348	sprintf (vFilename, "%s.tree%d.vstat", sumtParams.sumtOutfile, treeNo+1);
4349	sprintf (cFilename, "%s.tree%d.con.tre", sumtParams.sumtOutfile, treeNo+1);
4350	sprintf (tFilename, "%s.tree%d.trprobs", sumtParams.sumtOutfile, treeNo+1);
4351	}
4352	else
4353	{
4354	sprintf (pFilename, "%s.parts", sumtParams.sumtOutfile);
4355	sprintf (sFilename, "%s.tstat", sumtParams.sumtOutfile);
4356	sprintf (vFilename, "%s.vstat", sumtParams.sumtOutfile);
4357	sprintf (cFilename, "%s.con.tre", sumtParams.sumtOutfile);
4358	sprintf (tFilename, "%s.trprobs", sumtParams.sumtOutfile);
4359	}
4360
4361	/* open files checking for over-write as appropriate */
4362	if ((fpParts = OpenNewMBPrintFile(pFilename)) == NULL)
4363	return ERROR;
4364	if ((fpTstat = OpenNewMBPrintFile(sFilename)) == NULL)
4365	{
4366	SafeFclose (&fpParts);
4367	return ERROR;
4368	}
4369	if ((fpVstat = OpenNewMBPrintFile(vFilename)) == NULL)
4370	{
4371	SafeFclose (&fpParts);
4372	SafeFclose (&fpTstat);
4373	return ERROR;
4374	}
4375	if ((fpCon = OpenNewMBPrintFile(cFilename)) == NULL)
4376	{
4377	SafeFclose (&fpParts);
4378	SafeFclose (&fpTstat);
4379	SafeFclose (&fpVstat);
4380	return ERROR;
4381	}
4382	if (sumtParams.calcTreeprobs == YES)
4383	{
4384	if ((fpTrees = OpenNewMBPrintFile(tFilename)) == NULL)
4385	{
4386	SafeFclose (&fpParts);
4387	SafeFclose (&fpTstat);
4388	SafeFclose (&fpVstat);
4389	SafeFclose (&fpCon);
4390	return ERROR;
4391	}
4392	}
4393
4394	/* print #NEXUS if appropriate */
4395	fprintf (fpCon, "#NEXUS\n\n");
4396	if (sumtParams.calcTreeprobs == YES)
4397	fprintf (fpTrees, "#NEXUS\n\n");
4398
4399	/* print unique identifiers to each file */
4400	len = (int) strlen (stamp);
4401	if (len > 1)
4402	{
4403	fprintf (fpParts, "[ID: %s]\n", stamp);
4404	fprintf (fpTstat, "[ID: %s]\n", stamp);
4405	fprintf (fpVstat, "[ID: %s]\n", stamp);
4406	fprintf (fpCon, "[ID: %s]\n", stamp);
4407	if (sumtParams.calcTreeprobs == YES)
4408	fprintf (fpTrees, "[ID: %s]\n", stamp);
4409	}
4410
4411	/* Reset noWarn and autoOverwrite */
4412	if (treeNo == sumtParams.numTrees - 1)
4413	{
4414	noWarn = oldNoWarn;
4415	autoOverwrite = oldAutoOverwrite;
4416	}
4417
4418	return (NO_ERROR);
4419	}
4420
4421
4422
4423
4424
4425	void PartCtrUppass (PartCtr r, PartCtr uppass, int index)
4426
4427	{
4428	if (r != NULL)
4429	{
4430	uppass[(*index)++] = r;
4431
4432	PartCtrUppass (r->left, uppass, index);
4433	PartCtrUppass (r->right, uppass, index);
4434	}
4435	}
4436
4437
4438
4439
4440
4441	/* PrintBrlensToFile: Print brlens to file */
4442	int PrintBrlensToFile (PartCtr **treeParts, int numTreeParts, int treeNo)
4443
4444	{
4445	int i, j, runNo, numBrlens;
4446	char filename[100];
4447	PartCtr *x;
4448	FILE *fp;
4449
4450	/* set file name */
4451	if (sumtParams.numTrees > 1)
4452	sprintf (filename, "%s.tree%d.brlens", sumtParams.sumtOutfile, treeNo+1);
4453	else
4454	sprintf (filename, "%s.brlens", sumtParams.sumtOutfile);
4455
4456	/* Open file checking for over-write as appropriate */
4457	if ((fp = OpenNewMBPrintFile(filename)) == NULL)
4458	return ERROR;
4459
4460	/* count number of brlens to print */
4461	for (i=0; i<numTreeParts; i++)
4462	{
4463	if (treeParts[i]->totCount < sumtParams.brlensFreqDisplay)
4464	break;
4465	}
4466	numBrlens = i;
4467
4468	/* print header */
4469	for (i=0; i<numBrlens; i++)
4470	{
4471	MrBayesPrintf (fp, "v[%d]", i+1);
4472	if (i==numBrlens-1)
4473	MrBayesPrintf (fp, "\n");
4474	else
4475	MrBayesPrintf (fp, "\t");
4476	}
4477
4478	/* print values */
4479	for (i=0; i<numBrlens; i++)
4480	{
4481	x = treeParts[numBrlens];
4482	for (runNo=0; runNo<sumtParams.numRuns; runNo++)
4483	{
4484	MrBayesPrintf (fp, "%s", MbPrintNum (x->length[runNo][0]));
4485	for (j=1; j<x->count[i]; j++)
4486	{
4487	MrBayesPrintf (fp, "\t%s", MbPrintNum (x->length[runNo][j]));
4488	}
4489	}
4490	MrBayesPrintf (fp, "\n");
4491	}
4492
4493	return NO_ERROR;
4494	}
4495
4496
4497
4498
4499
4500	/* PrintConTree: Print consensus tree in standard format readable by TreeView, Paup etc */
4501	void PrintConTree (FILE fp, PolyTree t)
4502	{
4503	MrBayesPrintf (fp, " [Note: This tree contains information on the topology, \n");
4504	MrBayesPrintf (fp, " branch lengths (if present), and the probability\n");
4505	MrBayesPrintf (fp, " of the partition indicated by the branch.]\n");
4506	if (!strcmp(sumtParams.sumtConType, "Halfcompat"))
4507	MrBayesPrintf (fp, " tree con_50_majrule = ");
4508	else
4509	MrBayesPrintf (fp, " tree con_all_compat = ");
4510	WriteConTree (t->root, fp, YES);
4511	MrBayesPrintf (fp, ";\n");
4512	if (sumtParams.brlensDef == YES)
4513	{
4514	MrBayesPrintf (fp, "\n");
4515	MrBayesPrintf (fp, " [Note: This tree contains information only on the topology\n");
4516	MrBayesPrintf (fp, " and branch lengths (median of the posterior probability density).]\n");
4517	if (!strcmp(sumtParams.sumtConType, "Halfcompat"))
4518	MrBayesPrintf (fp, " tree con_50_majrule = ");
4519	else
4520	MrBayesPrintf (fp, " tree con_all_compat = ");
4521	WriteConTree (t->root, fp, NO);
4522	MrBayesPrintf (fp, ";\n");
4523	}
4524	}
4525
4526
4527
4528
4529
4530	/* PrintFigTreeConTree: Print consensus tree in rich format for FigTree */
4531	void PrintFigTreeConTree (FILE fp, PolyTree t, PartCtr **treeParts)
4532	{
4533	if (!strcmp(sumtParams.sumtConType, "Halfcompat"))
4534	MrBayesPrintf (fp, " tree con_50_majrule = ");
4535	else
4536	MrBayesPrintf (fp, " tree con_all_compat = ");
4537	if (t->isRooted == YES)
4538	MrBayesPrintf (fp, "[&R] ");
4539	else
4540	MrBayesPrintf (fp, "[&U] ");
4541
4542	WriteFigTreeConTree (t->root, fp, treeParts);
4543	MrBayesPrintf (fp, ";\n");
4544	}
4545
4546
4547
4548
4549
4550	void PrintFigTreeNodeInfo (FILE fp, PartCtr x, MrBFlt length)
4551
4552	{
4553	int i, postProbPercent, postProbSdPercent;
4554	MrBFlt *support, mean, var, min, max;
4555	Stat theStats;
4556
4557	support = SafeCalloc (sumtParams.numRuns, sizeof(MrBFlt));
4558	for (i=0; i<sumtParams.numRuns; i++)
4559	{
4560	support[i] = (MrBFlt) x->count[i] / (MrBFlt) sumtParams.numFileTreesSampled[i];
4561	}
4562	if (sumtParams.numRuns > 1)
4563	{
4564	MeanVariance (support, sumtParams.numRuns, &mean, &var);
4565	Range (support, sumtParams.numRuns, &min, &max);
4566	postProbPercent = (int) (100.0*mean + 0.5);
4567	postProbSdPercent = (int) (100.0 * sqrt(var) + 0.5);
4568	fprintf (fp, "[&prob=%.15le,prob_stddev=%.15le,prob_range={%.15le,%.15le},prob(percent)=\"%d\",prob+-sd=\"%d+-%d\"",
4569	mean, sqrt(var), min, max, postProbPercent, postProbPercent, postProbSdPercent);
4570	}
4571	else
4572	{
4573	postProbPercent = (int) (100.0*support[0] + 0.5);
4574	fprintf (fp, "[&prob=%.15le,prob(percent)=\"%d\"", support[0], postProbPercent);
4575	}
4576	if (sumtParams.isClock == YES)
4577	{
4578	GetSummary (x->height, sumtParams.numRuns, x->count, &theStats, sumtParams.HPD);
4579	if (sumtParams.HPD == YES)
4580	fprintf (fp, ",height_mean=%.15le,height_median=%.15le,height_95%%HPD={%.15le,%.15le}", theStats.mean, theStats.median, theStats.lower, theStats.upper);
4581	else
4582	fprintf (fp, ",height_mean=%.15le,height_median=%.15le,height_95%%CredInt={%.15le,%.15le}", theStats.mean, theStats.median, theStats.lower, theStats.upper);
4583	}
4584	if (sumtParams.isCalibrated == YES)
4585	{
4586	GetSummary (x->age, sumtParams.numRuns, x->count, &theStats, sumtParams.HPD);
4587	if (sumtParams.HPD == YES)
4588	fprintf (fp, ",age_mean=%.15le,age_median=%.15le,age_95%%HPD={%.15le,%.15le}", theStats.mean, theStats.median, theStats.lower, theStats.upper);
4589	else
4590	fprintf (fp, ",age_mean=%.15le,age_median=%.15le,age_95%%CredInt={%.15le,%.15le}", theStats.mean, theStats.median, theStats.lower, theStats.upper);
4591	}
4592	fprintf (fp, "]");
4593	if (length >= 0.0)
4594	fprintf (fp, ":%s", MbPrintNum(length));
4595	if (sumtParams.brlensDef == YES)
4596	{
4597	GetSummary (x->length, sumtParams.numRuns, x->count, &theStats, sumtParams.HPD);
4598	if (sumtParams.HPD == YES)
4599	fprintf (fp, "[&length_mean=%.15le,length_median=%.15le,length_95%%HPD={%.15le,%.15le}", theStats.mean, theStats.median, theStats.lower, theStats.upper);
4600	else
4601	fprintf (fp, "[&length_mean=%.15le,length_median=%.15le,length_95%%CredInt={%.15le,%.15le}", theStats.mean, theStats.median, theStats.lower, theStats.upper);
4602	}
4603	if (sumtParams.isClock == YES && sumtParams.isRelaxed == YES)
4604	{
4605	for (i=0; i<sumtParams.nBSets; i++)
4606	{
4607	GetSummary (x->bLen[i], sumtParams.numRuns, x->count, &theStats, sumtParams.HPD);
4608	if (sumtParams.HPD == YES)
4609	fprintf (fp, ",effectivebrlen%s_mean=%lf,effectivebrlen%s_median=%lf,effectivebrlen%s_95%%HPD={%lf,%lf}",
4610	sumtParams.tree->bSetName[i], theStats.mean,
4611	sumtParams.tree->bSetName[i], theStats.median,
4612	sumtParams.tree->bSetName[i], theStats.lower,
4613	theStats.upper);
4614	else
4615	fprintf (fp, ",effectivebrlen%s_mean=%lf,effectivebrlen%s_median=%lf,effectivebrlen%s_95%%CredInt={%lf,%lf}",
4616	sumtParams.tree->bSetName[i], theStats.mean,
4617	sumtParams.tree->bSetName[i], theStats.median,
4618	sumtParams.tree->bSetName[i], theStats.lower,
4619	theStats.upper);
4620	GetSummary (x->bRate[i], sumtParams.numRuns, x->count, &theStats, sumtParams.HPD);
4621	if (sumtParams.HPD == YES)
4622	fprintf (fp, ",rate%s_mean=%lf,rate%s_median=%lf,rate%s_95%%HPD={%lf,%lf}",
4623	sumtParams.tree->bSetName[i], theStats.mean,
4624	sumtParams.tree->bSetName[i], theStats.median,
4625	sumtParams.tree->bSetName[i], theStats.lower,
4626	theStats.upper);
4627	else
4628	fprintf (fp, ",rate%s_mean=%lf,rate%s_median=%lf,rate%s_95%%CredInt={%lf,%lf}",
4629	sumtParams.tree->bSetName[i], theStats.mean,
4630	sumtParams.tree->bSetName[i], theStats.median,
4631	sumtParams.tree->bSetName[i], theStats.lower,
4632	theStats.upper);
4633	}
4634	for (i=0; i<sumtParams.nESets; i++)
4635	{
4636	GetIntSummary (x->nEvents[i], sumtParams.numRuns, x->count, &theStats, sumtParams.HPD);
4637	if (sumtParams.HPD == YES)
4638	fprintf (fp, ",nEvents%s_mean=%lf,nEvents%s_median=%lf,nEvents%s_95%%HPD={%lf,%lf}",
4639	sumtParams.tree->eSetName[i], theStats.mean,
4640	sumtParams.tree->eSetName[i], theStats.median,
4641	sumtParams.tree->eSetName[i], theStats.lower,
4642	theStats.upper);
4643	else
4644	fprintf (fp, ",nEvents%s_mean=%lf,nEvents%s_median=%lf,nEvents%s_95%%CredInt={%lf,%lf}",
4645	sumtParams.tree->eSetName[i], theStats.mean,
4646	sumtParams.tree->eSetName[i], theStats.median,
4647	sumtParams.tree->eSetName[i], theStats.lower,
4648	theStats.upper);
4649	}
4650	}
4651	if (sumtParams.brlensDef == YES)
4652	fprintf (fp, "]");
4653
4654	free (support);
4655	}
4656
4657
4658
4659
4660
4661	/* PrintSumtTaxaInfo: Print information on pruned and absent taxa */
4662	void PrintSumtTaxaInfo (void)
4663	{
4664	int i, j, lineWidth, numExcludedTaxa, len;
4665	char tempStr[100];
4666
4667	/* print out information on absent taxa */
4668	numExcludedTaxa = 0;
4669	for (i=0; i<numTaxa; i++)
4670	if (sumtParams.absentTaxa[i] == YES)
4671	numExcludedTaxa++;
4672
4673	if (numExcludedTaxa > 0)
4674	{
4675	if (numExcludedTaxa == 1)
4676	MrBayesPrint ("%s The following taxon was absent from trees:\n", spacer);
4677	else
4678	MrBayesPrint ("%s The following %d taxa were absent from trees:\n", spacer, numExcludedTaxa);
4679	MrBayesPrint ("%s ", spacer);
4680	j = lineWidth = 0;
4681	for (i=0; i<numTaxa; i++)
4682	{
4683	if (sumtParams.absentTaxa[i] == YES)
4684	{
4685	j++;
4686	strcpy (tempStr, taxaNames[i]);
4687	len = (int) strlen(tempStr);
4688	lineWidth += len+2;
4689	if (lineWidth > 60)
4690	{
4691	MrBayesPrint ("\n%s ", spacer);
4692	lineWidth = 0;
4693	}
4694	if (numExcludedTaxa == 1)
4695	MrBayesPrint ("%s\n", tempStr);
4696	else if (numExcludedTaxa == 2 && j == 1)
4697	MrBayesPrint ("%s ", tempStr);
4698	else if (j == numExcludedTaxa)
4699	MrBayesPrint ("and %s\n", tempStr);
4700	else
4701	MrBayesPrint ("%s, ", tempStr);
4702	}
4703	}
4704	MrBayesPrint ("\n");
4705	}
4706
4707	/* print out information on pruned taxa */
4708	numExcludedTaxa = 0;
4709	for (i=0; i<numTaxa; i++)
4710	if (taxaInfo[i].isDeleted == YES && sumtParams.absentTaxa[i] == NO)
4711	numExcludedTaxa++;
4712
4713	if (numExcludedTaxa > 0)
4714	{
4715	if (numExcludedTaxa == 1)
4716	MrBayesPrint ("%s The following taxon was pruned from trees:\n", spacer);
4717	else
4718	MrBayesPrint ("%s The following %d taxa were pruned from trees:\n", spacer, numExcludedTaxa);
4719	MrBayesPrint ("%s ", spacer);
4720	j = lineWidth = 0;
4721	for (i=0; i<numTaxa; i++)
4722	{
4723	if (taxaInfo[i].isDeleted == YES && sumtParams.absentTaxa[i] == NO)
4724	{
4725	j++;
4726	strcpy (tempStr, taxaNames[i]);
4727	len = (int) strlen(tempStr);
4728	lineWidth += len+2;
4729	if (lineWidth > 60)
4730	{
4731	MrBayesPrint ("\n%s ", spacer);
4732	lineWidth = 0;
4733	}
4734	if (numExcludedTaxa == 1)
4735	MrBayesPrint ("%s\n", tempStr);
4736	else if (numExcludedTaxa == 2 && j == 1)
4737	MrBayesPrint ("%s ", tempStr);
4738	else if (j == numExcludedTaxa)
4739	MrBayesPrint ("and %s\n", tempStr);
4740	else
4741	MrBayesPrint ("%s, ", tempStr);
4742	}
4743	}
4744	MrBayesPrint ("\n");
4745	}
4746	}
4747
4748
4749
4750
4751
4752	/* Range: Determine range for a vector of MrBFlt values */
4753	void Range (MrBFlt vals, int nVals, MrBFlt min, MrBFlt *max)
4754
4755	{
4756	SortMrBFlt (vals, 0, nVals-1);
4757
4758	*min = vals[0];
4759	*max = vals[nVals-1];
4760
4761	}
4762
4763
4764
4765
4766
4767	/* ResetTaxonSet: Reset included taxa and local outgroup number */
4768	void ResetTaxonSet (void)
4769	{
4770	int i, j;
4771
4772	/* reset numLocalTaxa and localOutGroup */
4773	localOutGroup = 0;
4774	numLocalTaxa = 0;
4775	for (i=j=0; i<numTaxa; i++)
4776	{
4777	if (taxaInfo[i].isDeleted == NO)
4778	{
4779	if (i == outGroupNum)
4780	localOutGroup = numLocalTaxa;
4781	numLocalTaxa++;
4782	}
4783	}
4784
4785	}
4786
4787
4788
4789
4790
4791	void ResetTranslateTable (void)
4792	{
4793	int i;
4794
4795	for (i=0; i<numTranslates; i++)
4796	{
4797	free (transFrom[i]);
4798	free (transTo[i]);
4799	}
4800	free (transFrom);
4801	free (transTo);
4802	transFrom = NULL;
4803	transTo = NULL;
4804	numTranslates = 0;
4805	isTranslateDef = NO;
4806	isTranslateDiff = NO;
4807	}
4808
4809
4810
4811
4812
4813	int ShowConPhylogram (FILE fp, PolyTree t, int screenWidth)
4814
4815	{
4816
4817	int i, j, k, nLines, from, to, treeWidth=0, barLength, printExponential,
4818	precision, width, newPos, curPos, nTimes, numSpaces, maxLabelLength;
4819	char printLine, markLine, temp[30], *label;
4820	MrBFlt scale, f, scaleBar;
4821	PolyNode p, q;
4822
4823	/* set max label length */
4824	maxLabelLength = 20;
4825
4826	/* allocate space for label, printLine and markLine */
4827	printLine = (char ) SafeCalloc ((2screenWidth+2),sizeof(char));
4828	label = (char *) SafeCalloc (maxLabelLength+1, sizeof(char));
4829	if (!printLine \|\| !label)
4830	return ERROR;
4831	markLine = printLine + screenWidth + 1;
4832
4833	/* calculate scale */
4834	scale = 0.0;
4835	t->root->f = 0.0;
4836	for (i=t->nNodes-2; i>=0; i--)
4837	{
4838	p = t->allDownPass[i];
4839	/* find distance to root in relevant units */
4840	if (sumtParams.isClock == YES && sumtParams.isCalibrated == NO)
4841	p->f = t->root->depth - p->depth;
4842	else if (sumtParams.isClock == YES && sumtParams.isCalibrated == YES)
4843	p->f = t->root->age - p->age;
4844	else
4845	p->f = p->anc->f + p->length;
4846	if (p->left == NULL)
4847	{
4848	f = p->f / (screenWidth - Label(p,YES,NULL,maxLabelLength) - 2);
4849	if (f > scale)
4850	{
4851	scale = f;
4852	treeWidth = screenWidth - Label(p,YES,NULL,maxLabelLength) - 2;
4853	}
4854	}
4855	}
4856
4857	/* calculate x coordinates */
4858	for (i=0; i<t->nNodes; i++)
4859	{
4860	p = t->allDownPass[i];
4861	p->x = (int) (0.5 + (p->f / scale));
4862	}
4863
4864	/* calculate y coordinates and lines to print */
4865	for (i=nLines=0; i<t->nNodes; i++)
4866	{
4867	p = t->allDownPass[i];
4868	if (p->left != NULL)
4869	{
4870	/* internal node */
4871	for (q=p->left->sib; q->sib!=NULL; q=q->sib)
4872	;
4873	p->y = (int) (0.5 + ((p->left->y + q->y) / 2.0));
4874	}
4875	else
4876	{
4877	/* terminal node */
4878	p->y = nLines;
4879	nLines += 2;
4880	}
4881	}
4882
4883	/* print tree line by line */
4884
4885	for (i=0; i<nLines; i++)
4886	{
4887	MrBayesPrint ("%s ", spacer);
4888	/* empty printLine */
4889	for (j=0; j<screenWidth; j++)
4890	{
4891	printLine[j] = ' ';
4892	}
4893	printLine[j]='\0';
4894
4895	for (j=0; j<t->nNodes; j++)
4896	{
4897	p = t->allDownPass[j];
4898	if (p->y != i)
4899	continue;
4900
4901	/* this branch should be printed */
4902	/* add branch */
4903	if (p->anc == NULL)
4904	{
4905	/* this is the root of the whole tree */
4906	printLine[p->x] = '+';
4907	}
4908	else
4909	{
4910	/* this is an ordinary branch */
4911	to = p->x;
4912	from = p->anc->x;
4913	for (k=from+1; k<=to; k++)
4914	printLine[k] = '-';
4915	if (p == p->anc->left)
4916	{
4917	if (markLine[from] == 0)
4918	printLine[from] = '/';
4919	else
4920	printLine[from] = '\|';
4921	markLine[from] ++;
4922	}
4923	else if (p->sib == NULL)
4924	{
4925	if (markLine[from] == 1)
4926	printLine[from] = '\\';
4927	else
4928	printLine[from] = '\|';
4929	markLine[from] --;
4930	}
4931	if (p->left!=NULL)
4932	{
4933	if (from != to)
4934	printLine[to] = '+';
4935	else
4936	printLine[to] = '\|';
4937	}
4938	else
4939	{
4940	/* add label if the branch is terminal */
4941	Label(p,YES,label,maxLabelLength);
4942	sprintf(printLine+to+2,"%s", label);
4943	}
4944	}
4945	}
4946
4947	/* check for cross branches */
4948	for (j=0; j<screenWidth; j++)
4949	{
4950	if (markLine[j] >= 1 && printLine[j] == ' ')
4951	printLine[j] = '\|';
4952	}
4953	MrBayesPrintf (fp, "%s\n",printLine);
4954	}
4955
4956	/* print scale */
4957	k = (int) (floor (log10 (scale * 80)));
4958	scaleBar = pow (10, k);
4959	barLength = (int) (scaleBar / scale);
4960	if (barLength > 80)
4961	{
4962	barLength /= 10;
4963	scaleBar /= 10.0;
4964	}
4965	else if (barLength > 40)
4966	{
4967	barLength /= 5;
4968	scaleBar /= 5.0;
4969	}
4970	else if (barLength > 16)
4971	{
4972	barLength /= 2;
4973	scaleBar /= 2.0;
4974	}
4975
4976	if (t->isClock == YES)
4977	{
4978	MrBayesPrint ("%s ", spacer);
4979	for (i=0; i<treeWidth; i++)
4980	printLine[i] = '-';
4981	nTimes = (int) (treeWidth / (scaleBar / scale));
4982	for (i=0; i<=nTimes; i++)
4983	printLine[treeWidth - (int)(i*scaleBar/scale)] = '\|';
4984	MrBayesPrint ("%s\n", printLine);
4985
4986	MrBayesPrint ("%s ", spacer);
4987	f = treeWidth * scale;
4988	if (f >= 1000.0 \|\| f < 0.10)
4989	{
4990	printExponential = YES;
4991	precision = 0;
4992	}
4993	else
4994	{
4995	printExponential = NO;
4996	precision = 2 - (int) (log10 (f));
4997	}
4998
4999	curPos = 0;
5000	f = nTimes * scaleBar;
5001	while (curPos < treeWidth)
5002	{
5003	/* print the number */
5004	if (printExponential == YES)
5005	sprintf (temp, "%.2e", f);
5006	else
5007	sprintf (temp, "%.*lf", precision, f);
5008
5009	/* room to print ? if so, print */
5010	width = (int) strlen (temp);
5011	newPos = treeWidth - (int) (nTimes * (scaleBar / scale));
5012	numSpaces = newPos - width / 2 - curPos;
5013	if (numSpaces >= 0 \|\| (numSpaces >= -2 && curPos == 0))
5014	{
5015	while (numSpaces > 0)
5016	{
5017	printLine[curPos++] = ' ';
5018	numSpaces--;
5019	}
5020	for (i=0; temp[i]!='\0'; i++)
5021	printLine[curPos++] = temp[i];
5022	}
5023
5024	/* get new number */
5025	f -= scaleBar;
5026	nTimes--;
5027	}
5028
5029	MrBayesPrint ("%s\n", printLine);
5030
5031	if (sumtParams.isCalibrated == YES)
5032	MrBayesPrint ("\n%s [User-defined time units]\n\n", spacer);
5033	else
5034	MrBayesPrint ("\n%s [Expected changes per site]\n\n", spacer);
5035	}
5036	else
5037	{
5038	MrBayesPrintf (fp, "%s \|", spacer);
5039	for (i=0; i<barLength-1; i++)
5040	MrBayesPrintf (fp, "-");
5041	MrBayesPrintf (fp, "\| %1.3lf expected changes per site\n\n", scaleBar);
5042	}
5043
5044	free (printLine);
5045
5046	return NO_ERROR;
5047
5048	}
5049
5050
5051
5052	int ShowConTree (FILE fp, PolyTree t, int screenWidth, int showSupport)
5053
5054	{
5055
5056	int i, j, k, treeWidth, minBranchLength, maxWidth, isTreeDivided,
5057	printWidth, nLines, nodesToBePrinted, from, to, maxLabelLength,
5058	maxLength;
5059	char printLine, markLine, temp[20], *label;
5060	PolyNode p=NULL, q;
5061
5062	maxLength = 20; /* max length of label */
5063	minBranchLength = 5; /* min length of branch in tree */
5064	isTreeDivided = NO;
5065
5066	/* allocate space for printLine, markLine and label */
5067	printLine = (char ) SafeCalloc (maxLength+1+(2screenWidth+2),sizeof(char));
5068	if (!printLine)
5069	return ERROR;
5070	markLine = printLine + screenWidth + 1;
5071	label = markLine + screenWidth + 1;
5072
5073	/* get fresh internal node indices */
5074	k = t->nNodes - t->nIntNodes;
5075	for (i=0; i<t->nIntNodes; i++)
5076	{
5077	p = t->intDownPass[i];
5078	p->index = k++;
5079	}
5080
5081	/* calculate max length of labels including taxon index number */
5082	maxLabelLength = 0;
5083	for (i=0; i<t->nNodes; i++)
5084	{
5085	p = t->allDownPass[i];
5086	if (p->left == NULL)
5087	{
5088	j = Label(p,YES,NULL,maxLength);
5089	if (j > maxLabelLength)
5090	maxLabelLength = j;
5091	}
5092	}
5093
5094	/* make sure label can hold an interior node index number */
5095	j = (int) (3.0 + log10((MrBFlt)t->nNodes));
5096	maxLabelLength = (maxLabelLength > j ? maxLabelLength : j);
5097
5098	/* calculate remaining screen width for tree
5099	and maxWidth in terms of branches */
5100	treeWidth = screenWidth - maxLabelLength - 1;
5101	maxWidth = treeWidth / minBranchLength;
5102
5103	/* unmark whole tree */
5104	for (i=0; i<t->nNodes; i++)
5105	t->allDownPass[i]->mark = 0;
5106	nodesToBePrinted = t->nNodes;
5107
5108	while (nodesToBePrinted > 0)
5109	{
5110	/* count depth of nodes in unprinted tree */
5111	for (i=0; i<t->nNodes; i++)
5112	{
5113	p = t->allDownPass[i];
5114	if (p->mark == 0) /* the node has not been printed yet */
5115	{
5116	p->x = 0;
5117	/* if it is an interior node in the tree that will be printed
5118	compute the depth of the node */
5119	if (p->left != NULL && p->left->mark == 0)
5120	{
5121	for (q = p->left; q!=NULL; q=q->sib)
5122	{
5123	if (q->x > p->x)
5124	p->x = q->x;
5125	}
5126	p->x++;
5127	/* break when root of print subtree has been found */
5128	if (p->x >= maxWidth)
5129	break;
5130	}
5131	}
5132	}
5133
5134	/* if internal node then find largest nonprinted subtree among descendant nodes */
5135	if (p->anc != NULL)
5136	{
5137	for (q=p->left; q!=NULL; q=q->sib)
5138	{
5139	if (q->x == p->x - 1 && q->mark == 0)
5140	p = q;
5141	}
5142	MrBayesPrintf (fp, "%s Subtree rooted at node %d:\n\n", spacer, p->index);
5143	isTreeDivided = YES;
5144	}
5145	else if (isTreeDivided == YES)
5146	MrBayesPrintf (fp, "%s Root part of tree:\n\n", spacer);
5147
5148	/* mark subtree for printing and
5149	translate x coordinates from depth to position */
5150	if (p->anc == NULL)
5151	printWidth = p->x;
5152	else
5153	printWidth = p->x + 1;
5154	p->mark = 1;
5155	p->x = (int) (treeWidth - 0.5 - ((treeWidth - 1) * (p->x / (MrBFlt) printWidth)));
5156	for (i=t->nNodes-2; i>=0; i--)
5157	{
5158	p = t->allDownPass[i];
5159	if (p->mark == 0 && p->anc->mark == 1)
5160	{
5161	p->mark = 1;
5162	p->x = (int) (treeWidth - 0.5 - ((treeWidth - 1) * (p->x / (MrBFlt) printWidth)));
5163	}
5164	}
5165
5166	/* calculate y coordinates of nodes to be printed and lines to print */
5167	for (i=nLines=0; i<t->nNodes; i++)
5168	{
5169	p = t->allDownPass[i];
5170	if (p->mark == 1)
5171	{
5172	if (p->left != NULL && p->left->mark == 1)
5173	{
5174	/* internal node */
5175	for (q=p->left->sib; q->sib!=NULL; q=q->sib)
5176	;
5177	p->y = (int) (0.5 + ((p->left->y + q->y) / 2.0));
5178	}
5179	else
5180	{
5181	/* terminal node */
5182	p->y = nLines;
5183	nLines += 2;
5184	}
5185	}
5186	}
5187
5188	/* print subtree line by line */
5189	for (i=0; i<nLines; i++)
5190	{
5191	MrBayesPrintf (fp, "%s ", spacer);
5192	/* empty printLine */
5193	for (j=0; j<screenWidth; j++)
5194	{
5195	printLine[j] = ' ';
5196	}
5197	printLine[j]='\0';
5198
5199	for (j=0; j<t->nNodes; j++)
5200	{
5201	p = t->allDownPass[j];
5202	if (p->mark != 1 \|\| p->y != i)
5203	continue;
5204
5205	/* this branch should be printed
5206	add label if the branch is terminal in tree to be printed */
5207	if (p->left == NULL)
5208	{
5209	Label (p,YES,label,maxLength);
5210	sprintf(printLine+treeWidth+1,"%s", label);
5211	}
5212	else if (p->left->mark == 2)
5213	sprintf(printLine+treeWidth+1,"(%d)", p->index);
5214
5215	/* add branch */
5216	if (p->anc == NULL)
5217	{
5218	/* this is the root of the whole tree */
5219	printLine[p->x] = '+';
5220	nodesToBePrinted--;
5221	}
5222	else if (p->anc->mark == 0)
5223	{
5224	/* this is a root of a subtree
5225	this branch will have to be printed again so do
5226	not decrease nodesToBePrinted */
5227	to = p->x;
5228	from = 0;
5229	for (k=from; k<to; k++)
5230	printLine[k] = '-';
5231	printLine[to] = '+';
5232	if (showSupport == YES)
5233	sprintf(temp, "%d", (int) (p->support*100.0 + 0.5));
5234	else
5235	*temp='\0';
5236	from = (int)(from + 1.5 + ((to - from - 1 - strlen(temp)) / 2.0));
5237	for (k=0; temp[k]!='\0'; k++)
5238	printLine[from++] = temp[k];
5239	}
5240	else
5241	{
5242	/* this is an ordinary branch */
5243	to = p->x;
5244	from = p->anc->x;
5245	for (k=from+1; k<=to; k++)
5246	printLine[k] = '-';
5247	if (p == p->anc->left)
5248	{
5249	printLine[from] = '/';
5250	markLine[from] = 1;
5251	}
5252	else if (p->sib == NULL)
5253	{
5254	printLine[from] = '\\';
5255	markLine[from] = 0;
5256	}
5257	if (p->left!=NULL && p->left->mark!=2)
5258	{
5259	printLine[to] = '+';
5260	if (showSupport == YES)
5261	sprintf(temp, "%d", (int) (p->support*100.0 + 0.5));
5262	else
5263	*temp='\0';
5264	from = (int)(from + 1.5 + ((to - from - 1 - strlen(temp)) / 2.0));
5265	for (k=0; temp[k]!='\0'; k++)
5266	printLine[from++] = temp[k];
5267	}
5268	nodesToBePrinted--;
5269	}
5270	}
5271
5272	/* check for cross branches */
5273	for (j=0; j<treeWidth; j++)
5274	{
5275	if (markLine[j] == 1 && printLine[j] == ' ')
5276	printLine[j] = '\|';
5277	}
5278
5279	MrBayesPrintf(fp, "%s\n",printLine);
5280	}
5281
5282	/* mark printed branches */
5283	for (i=0; i<t->nNodes; i++)
5284	{
5285	p = t->allDownPass[i];
5286	if (p->mark == 1)
5287	{
5288	if (p->anc == NULL)
5289	p->mark = 2;
5290	else if (p->anc->mark == 0)
5291	p->mark = 0; /* this branch will have to be printed again */
5292	else
5293	p->mark = 2;
5294	}
5295	}
5296
5297	} /* next subtree */
5298
5299	free (printLine);
5300
5301	return NO_ERROR;
5302
5303	}
5304
5305
5306
5307
5308
5309	void ShowParts (FILE fp, SafeLong p, int nTaxaToShow)
5310
5311	{
5312
5313	int i;
5314	SafeLong x, y;
5315
5316	for (i=0; i<nTaxaToShow; i++)
5317	{
5318	y = p[i / nBitsInALong];
5319	x = 1 << (i % nBitsInALong);
5320	if ((x & y) == 0)
5321	MrBayesPrintf (fp, ".");
5322	else
5323	MrBayesPrintf (fp, "*");
5324	}
5325
5326	}
5327
5328
5329
5330
5331
5332	void ShowSomeParts (FILE fp, SafeLong p, int offset, int nTaxaToShow)
5333
5334	{
5335
5336	int i;
5337	SafeLong x, y;
5338
5339	for (i=offset; i<offset+nTaxaToShow; i++)
5340	{
5341	y = p[i / nBitsInALong];
5342	x = 1 << (i % nBitsInALong);
5343	if ((x & y) == 0)
5344	MrBayesPrintf (fp, ".");
5345	else
5346	MrBayesPrintf (fp, "*");
5347	}
5348	}
5349
5350
5351
5352
5353
5354	void SortPartCtr (PartCtr **item, int left, int right)
5355
5356	{
5357
5358	register int i, j;
5359	PartCtr *tempPartCtr;
5360	int x;
5361
5362	assert (left >= 0);
5363	assert (right >= 0);
5364
5365	i = left;
5366	j = right;
5367	x = item[(left+right)/2]->totCount;
5368	do
5369	{
5370	while (item[i]->totCount > x && i < right)
5371	i++;
5372	while (x > item[j]->totCount && j > left)
5373	j--;
5374	if (i <= j)
5375	{
5376	tempPartCtr = item[i];
5377	item[i] = item[j];
5378	item[j] = tempPartCtr;
5379
5380	i++;
5381	j--;
5382	}
5383	} while (i <= j);
5384	if (left < j)
5385	SortPartCtr (item, left, j);
5386	if (i < right)
5387	SortPartCtr (item, i, right);
5388	}
5389
5390
5391
5392
5393
5394	void SortTerminalPartCtr (PartCtr **item, int len)
5395
5396	{
5397
5398	register int i, j, maxCount;
5399	PartCtr *temp;
5400
5401	maxCount = item[0]->totCount;
5402
5403	/* put root first */
5404	for (i=0; item[i]->totCount == maxCount; i++)
5405	if (NumBits(item[i]->partition, sumtParams.SafeLongsNeeded) == sumtParams.numTaxa)
5406	break;
5407
5408	if (i!=0)
5409	{
5410	temp = item[0];
5411	item[0] = item[i];
5412	item[i] = temp;
5413	}
5414
5415	/* then find terminals in index order */
5416	for (i=1; i<=sumtParams.numTaxa; i++)
5417	{
5418	for (j=i; item[j]->totCount == maxCount && j<len; j++)
5419	if (NumBits(item[j]->partition, sumtParams.SafeLongsNeeded) == 1 &&
5420	FirstTaxonInPartition(item[j]->partition, sumtParams.SafeLongsNeeded) == i-1)
5421	break;
5422
5423	if (j!=i)
5424	{
5425	temp = item[i];
5426	item[i] = item[j];
5427	item[j] = temp;
5428	}
5429	}
5430	}
5431
5432
5433
5434
5435
5436	void SortTreeCtr (TreeCtr **item, int left, int right)
5437
5438	{
5439
5440	register int i, j;
5441	TreeCtr *tempTreeCtr;
5442	int x;
5443
5444	i = left;
5445	j = right;
5446	x = item[(left+right)/2]->count;
5447	do
5448	{
5449	while (item[i]->count > x && i < right)
5450	i++;
5451	while (x > item[j]->count && j > left)
5452	j--;
5453	if (i <= j)
5454	{
5455	tempTreeCtr = item[i];
5456	item[i] = item[j];
5457	item[j] = tempTreeCtr;
5458
5459	i++;
5460	j--;
5461	}
5462	} while (i <= j);
5463	if (left < j)
5464	SortTreeCtr (item, left, j);
5465	if (i < right)
5466	SortTreeCtr (item, i, right);
5467	}
5468
5469
5470
5471
5472
5473	/* StoreSumtTree: Store tree in treeList in packed format */
5474	int StoreSumtTree (PackedTree treeList, int index, PolyTree t)
5475	{
5476	int orderLen, numBrlens;
5477
5478	assert(treeList[index].brlens == NULL);
5479	assert(treeList[index].order == NULL);
5480
5481	/* get tree dimensions */
5482	numBrlens = t->nNodes - 1;
5483	orderLen = t->nIntNodes - 1;
5484
5485	/* allocate space */
5486	treeList[index].brlens = (MrBFlt *) SafeCalloc (numBrlens, sizeof(MrBFlt));
5487	treeList[index].order = (int *) SafeCalloc (orderLen, sizeof(MrBFlt));
5488	if (!treeList[index].order \|\| !treeList[index].brlens)
5489	{
5490	MrBayesPrint ("%s Could not store packed representation of tree '%s'\n", spacer, t->name);
5491	return (ERROR);
5492	}
5493
5494	/* store tree */
5495	if (t->isRooted == YES)
5496	StoreRPolyTree (t, treeList[index].order, treeList[index].brlens);
5497	else
5498	StoreUPolyTree (t, treeList[index].order, treeList[index].brlens);
5499
5500	return (NO_ERROR);
5501	}
5502
5503
5504
5505
5506
5507	/* TreeCtrUppass: extract TreeCtr nodes in uppass sequence */
5508	void TreeCtrUppass (TreeCtr r, TreeCtr uppass, int index)
5509
5510	{
5511	if (r != NULL)
5512	{
5513	uppass[(*index)++] = r;
5514
5515	TreeCtrUppass (r->left, uppass, index);
5516	TreeCtrUppass (r->right, uppass, index);
5517	}
5518	}
5519
5520
5521
5522
5523
5524	int TreeProb (void)
5525
5526	{
5527
5528	int i, j, num, nInSets[5];
5529	MrBFlt treeProb, cumTreeProb;
5530	TreeCtr **trees;
5531	Tree *theTree;
5532
5533	/* check if we need to do this */
5534	if (sumtParams.calcTreeprobs == NO)
5535	return (NO_ERROR);
5536
5537	MrBayesPrint ("%s Calculating tree probabilities...\n\n", spacer);
5538
5539	/* allocate space for tree counters and trees */
5540	trees = (TreeCtr *) SafeCalloc ((size_t)numUniqueTreesFound, sizeof(TreeCtr ));
5541	theTree = AllocateTree (sumtParams.numTaxa);
5542	if (!trees \|\| !theTree)
5543	{
5544	MrBayesPrint ("%s Problem allocating trees or theTree in TreeProb\n", spacer);
5545	return (ERROR);
5546	}
5547
5548	/* extract trees */
5549	i = 0;
5550	TreeCtrUppass (treeCtrRoot, trees, &i);
5551
5552	/* sort trees */
5553	SortTreeCtr (trees, 0, numUniqueTreesFound-1);
5554
5555	/* set basic params in receiving tree */
5556	theTree->isRooted = sumtParams.isRooted;
5557	if (theTree->isRooted)
5558	{
5559	theTree->nNodes = 2 * sumtParams.numTaxa;
5560	theTree->nIntNodes = sumtParams.numTaxa - 1;
5561	}
5562	else
5563	{
5564	theTree->nNodes = 2 * sumtParams.numTaxa - 2;
5565	theTree->nIntNodes = sumtParams.numTaxa - 2;
5566	}
5567
5568	/* show tree data */
5569	cumTreeProb = 0.0;
5570	nInSets[0] = nInSets[1] = nInSets[2] = nInSets[3] = nInSets[4] = 0;
5571	for (num=0; num<numUniqueTreesFound; num++) /* loop over all of the trees that were found */
5572	{
5573	/* get probability of tree */
5574	treeProb = (MrBFlt)trees[num]->count / (MrBFlt)sumtParams.numTreesSampled;
5575	cumTreeProb += treeProb;
5576	if (cumTreeProb >= 0.0 && cumTreeProb < 0.5)
5577	nInSets[0]++;
5578	else if (cumTreeProb >= 0.5 && cumTreeProb < 0.9)
5579	nInSets[1]++;
5580	else if (cumTreeProb >= 0.9 && cumTreeProb < 0.95)
5581	nInSets[2]++;
5582	else if (cumTreeProb >= 0.95 && cumTreeProb < 0.99)
5583	nInSets[3]++;
5584	else
5585	nInSets[4]++;
5586
5587	/* draw tree to stdout */
5588	if (theTree->isRooted == YES)
5589	RetrieveRTopology (theTree, trees[num]->order);
5590	else
5591	RetrieveUTopology (theTree, trees[num]->order);
5592	if (sumtParams.showSumtTrees == YES)
5593	{
5594	MrBayesPrint ("\n%s Tree %d (p = %1.3lf, P = %1.3lf):\n\n", spacer, num+1, treeProb, cumTreeProb);
5595	ShowTree (theTree);
5596	}
5597
5598	/* draw tree to file */
5599	if (num == 0)
5600	{
5601	MrBayesPrintf (fpTrees, "[This file contains the trees that were found during the MCMC\n");
5602	MrBayesPrintf (fpTrees, "search, sorted by posterior probability. \"p\" indicates the\n");
5603	MrBayesPrintf (fpTrees, "posterior probability of the tree whereas \"P\" indicates the\n");
5604	MrBayesPrintf (fpTrees, "cumulative posterior probability.]\n\n");
5605	MrBayesPrintf (fpTrees, "begin trees;\n");
5606	MrBayesPrintf (fpTrees, " translate\n");
5607	j = 0;
5608	for (i=0; i<numTaxa; i++)
5609	{
5610	if (taxaInfo[i].isDeleted == NO && sumtParams.absentTaxa[i] == NO)
5611	{
5612	if (j+1 == sumtParams.numTaxa)
5613	MrBayesPrintf (fpTrees, " %2d %s;\n", j+1, taxaNames[i]);
5614	else
5615	MrBayesPrintf (fpTrees, " %2d %s,\n", j+1, taxaNames[i]);
5616	j++;
5617	}
5618	}
5619	}
5620	MrBayesPrintf (fpTrees, " tree tree_%d [p = %1.3lf, P = %1.3lf] = [&W %1.6lf] ", num+1, treeProb, cumTreeProb, treeProb);
5621	WriteTopologyToFile (fpTrees, theTree->root->left, theTree->isRooted);
5622	MrBayesPrintf (fpTrees, ";\n");
5623	if (num == numUniqueTreesFound - 1)
5624	MrBayesPrintf (fpTrees, "end;\n");
5625	}
5626
5627	/* print out general information on credible sets of trees */
5628	i = nInSets[0] + nInSets[1] + nInSets[2] + nInSets[3] + nInSets[4];
5629	MrBayesPrint ("%s Credible sets of trees (%d tree%s sampled):\n", spacer, i, i > 1 ? "s" : "");
5630	i = nInSets[0] + 1;
5631	if (i > 1)
5632	MrBayesPrint ("%s 50 %% credible set contains %d trees\n", spacer, i);
5633	i += nInSets[1];
5634	if (i > 1)
5635	MrBayesPrint ("%s 90 %% credible set contains %d trees\n", spacer, i);
5636	i += nInSets[2];
5637	if (i > 1)
5638	MrBayesPrint ("%s 95 %% credible set contains %d trees\n", spacer, i);
5639	i += nInSets[3];
5640	MrBayesPrint ("%s 99 %% credible set contains %d tree%s\n\n", spacer, i, i > 1 ? "s" : "");
5641
5642	/* free memory */
5643	free (trees);
5644
5645	return (NO_ERROR);
5646	}
5647
5648
5649
5650
5651
5652	void WriteConTree (PolyNode p, FILE fp, int showSupport)
5653
5654	{
5655
5656	PolyNode *q;
5657
5658	if (p->anc != NULL)
5659	if (p->anc->left == p)
5660	fprintf (fp, "(");
5661
5662	for (q = p->left; q != NULL; q = q->sib)
5663	{
5664	if (q->anc->left != q) /* Note that q->anc always exists (it is p) */
5665	fprintf (fp, ",");
5666	WriteConTree (q, fp, showSupport);
5667	}
5668	if (p->left == NULL)
5669	{
5670	if (sumtParams.brlensDef == YES)
5671	{
5672	if (sumtParams.isClock == NO)
5673	fprintf (fp, "%d:%s", p->index+1, MbPrintNum(p->length));
5674	else
5675	fprintf (fp, "%d:%s", p->index+1, MbPrintNum(p->anc->depth - p->depth));
5676	}
5677	else
5678	fprintf (fp, "%d", p->index+1);
5679	}
5680
5681	if (p->sib == NULL && p->anc != NULL)
5682	{
5683	if (p->anc->anc != NULL)
5684	{
5685	if (sumtParams.brlensDef == YES && showSupport == NO)
5686	{
5687	if (sumtParams.isClock == NO)
5688	fprintf (fp, "):%s", MbPrintNum(p->anc->length));
5689	else
5690	fprintf (fp, "):%s", MbPrintNum(p->anc->anc->depth - p->anc->depth));
5691	}
5692	else if (sumtParams.brlensDef == NO && showSupport == YES)
5693	fprintf (fp, ")%1.3lf", p->anc->support);
5694	else if (sumtParams.brlensDef == YES && showSupport == YES)
5695	{
5696	if (sumtParams.isClock == NO)
5697	fprintf (fp, ")%1.3lf:%s", p->anc->support, MbPrintNum(p->anc->length));
5698	else
5699	fprintf (fp, ")%1.3lf:%s", p->anc->support, MbPrintNum(p->anc->anc->depth - p->anc->depth));
5700	}
5701	else
5702	fprintf (fp, ")");
5703	}
5704	else
5705	fprintf (fp, ")");
5706	}
5707	}
5708
5709
5710
5711
5712
5713	/* WriteFigTreeConTree: Include rich information for each node in a consensus tree */
5714	void WriteFigTreeConTree (PolyNode p, FILE fp, PartCtr **treeParts)
5715
5716	{
5717
5718	PolyNode *q;
5719
5720	if (p->left == NULL)
5721	{
5722	fprintf (fp, "%d", p->index+1);
5723	if (sumtParams.isClock == NO)
5724	PrintFigTreeNodeInfo(fp, treeParts[p->partitionIndex], p->length);
5725	else if (sumtParams.isCalibrated == YES)
5726	PrintFigTreeNodeInfo(fp, treeParts[p->partitionIndex], p->anc->age - p->age);
5727	else
5728	PrintFigTreeNodeInfo(fp, treeParts[p->partitionIndex], p->anc->depth - p->depth);
5729	}
5730	else
5731	{
5732	fprintf (fp, "(");
5733	for (q = p->left; q != NULL; q = q->sib)
5734	{
5735	WriteFigTreeConTree (q, fp, treeParts);
5736	if (q->sib != NULL)
5737	fprintf (fp, ",");
5738	}
5739	fprintf (fp, ")");
5740	if (p->partitionIndex >= 0 && p->partitionIndex < numUniqueSplitsFound)
5741	{
5742	if (p->anc == NULL)
5743	PrintFigTreeNodeInfo(fp,treeParts[p->partitionIndex], -1.0);
5744	else if (sumtParams.isClock == NO)
5745	PrintFigTreeNodeInfo(fp, treeParts[p->partitionIndex], p->length);
5746	else if (sumtParams.isCalibrated == YES)
5747	PrintFigTreeNodeInfo(fp, treeParts[p->partitionIndex], p->anc->age - p->age);
5748	else
5749	PrintFigTreeNodeInfo(fp, treeParts[p->partitionIndex], p->anc->depth - p->depth);
5750	}
5751	}
5752	}
5753

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

Download in other formats: