7 * Dept. Integrative Biology
8 * University of California, Berkeley
9 * Berkeley, CA 94720-3140
13 * Swedish Museum of Natural History
15 * SE-10405 Stockholm, SWEDEN
16 * fredrik.ronquist@nrm.se
18 * With important contributions by
20 * Paul van der Mark (paulvdm@sc.fsu.edu)
21 * Maxim Teslenko (maxim.teslenko@nrm.se)
23 * and by many users (run 'acknowledgments' to see more info)
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.
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).
42 #if defined(__MWERKS__)
46 const char* const svnRevisionSumptC = "$Rev: 1072 $"; /* Revision keyword which is expended/updated by svn on each commit/update */
48 typedef struct partctr
50 struct partctr *left, *right;
57 int ***nEvents; /* nEvents[0,nESets][0,numRuns][0,count[RunID]] */
58 MrBFlt ***bRate; /* bRate [0,nBSets][0,numRuns][0,count[RunID]] */
59 MrBFlt ***bLen; /* bLen [0,nBSets][0,numRuns][0,count[RunID]] */
60 MrBFlt **popSize; /* popSize[0,numRuns][0,count[RunID]] */
64 typedef struct treectr
66 struct treectr *left, *right;
74 int longestLineLength;
76 int lastTreeBlockBegin;
78 int numTreesInLastBlock;
82 #define ALLOC_LEN 100 /* number of values to allocate each time in partition counter nodes */
84 #if defined (PRINT_RATEMUL_CPP)
85 FILE *rateMultfp=NULL;
90 extern int inSumtCommand;
91 extern int inComparetreeCommand;
92 extern int DoUserTree (void);
93 extern int DoUserTreeParm (char *parmName, char *tkn);
94 extern int SafeFclose(FILE **);
96 extern int SetUpPartitionCounters (void);
97 extern int AddTreeToPartitionCounters (Tree *tree, int treeId, int runId);
98 extern void CalcTopoConvDiagn2 (int *nTrees);
99 extern void FreeChainMemory (void);
101 /* local prototypes */
102 int CompareModelProbs (const void *x, const void *y);
103 int PrintModelStats (char *fileName, char **headerNames, int nHeaders, ParameterSample *parameterSamples, int nRuns, int nSamples);
104 int PrintOverlayPlot (MrBFlt **xVals, MrBFlt **yVals, int nRows, int startingFrom, int nSamples);
105 int PrintParamStats (char *fileName, char **headerNames, int nHeaders, ParameterSample *parameterSamples, int nRuns, int nSamples);
106 void PrintPlotHeader (void);
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, BitsLong *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);
140 void WriteConTree (PolyNode *p, FILE *fp, int showSupport);
141 void WriteFigTreeConTree (PolyNode *p, FILE *fp, PartCtr **treeParts);
143 /* local (to this file) */
144 static int numUniqueSplitsFound, numUniqueTreesFound, numPackedTrees[2], numAsterices; /* length of local to this file variables */
145 static FILE *fpParts=NULL, *fpTstat=NULL, *fpVstat, *fpCon=NULL, *fpTrees=NULL, *fpDists=NULL; /* file pointers */
146 static PartCtr *partCtrRoot = NULL; /* binary tree for holding splits info */
147 static TreeCtr *treeCtrRoot = NULL; /* binary tree for holding unique tree info */
148 static PackedTree *packedTreeList[2]; /* list of trees in packed format */
151 /* AllocateParameterSamples: Allocate space for parameter samples */
152 int AllocateParameterSamples (ParameterSample **parameterSamples, int numRuns, int numRows, int numColumns)
156 (*parameterSamples) = (ParameterSample *) SafeCalloc (numColumns, sizeof(ParameterSample));
157 if (!(*parameterSamples))
159 (*parameterSamples)[0].values = (MrBFlt **) SafeCalloc (numColumns * numRuns, sizeof (MrBFlt *));
160 if (!((*parameterSamples)[0].values))
162 FreeParameterSamples(*parameterSamples);
165 (*parameterSamples)[0].values[0] = (MrBFlt *) SafeCalloc (numColumns * numRuns * numRows, sizeof (MrBFlt));
166 for (i=1; i<numColumns; i++)
167 (*parameterSamples)[i].values = (*parameterSamples)[0].values + i*numRuns;
168 for (i=1; i<numRuns; i++)
169 (*parameterSamples)[0].values[i] = (*parameterSamples)[0].values[0] + i*numRows;
170 for (i=1; i<numColumns; i++)
172 for (j=0; j<numRuns; j++)
174 (*parameterSamples)[i].values[j] = (*parameterSamples)[0].values[0] + i*numRuns*numRows + j*numRows;
182 /** Compare function (ModelProb) for qsort. Note reverse sort order (from larger to smaller probs) */
183 int CompareModelProbs (const void *x, const void *y) {
185 if ((*((ModelProb *)(x))).prob > (*((ModelProb *)(y))).prob)
187 else if ((*((ModelProb *)(x))).prob < (*((ModelProb *)(y))).prob)
196 int i, n, nHeaders=0, numRows, numColumns, numRuns, whichIsX, whichIsY,
197 unreliable, oneUnreliable, burnin, longestHeader, len;
198 MrBFlt mean, harm_mean;
199 char **headerNames=NULL, temp[120];
200 SumpFileInfo fileInfo, firstFileInfo;
201 ParameterSample *parameterSamples=NULL;
204 # if defined (MPI_ENABLED)
209 /* tell user we are ready to go */
210 if (sumpParams.numRuns == 1)
211 MrBayesPrint ("%s Summarizing parameters in file %s.p\n", spacer, sumpParams.sumpFileName);
212 else if (sumpParams.numRuns == 2)
213 MrBayesPrint ("%s Summarizing parameters in files %s.run1.p and %s.run2.p\n", spacer, sumpParams.sumpFileName, sumpParams.sumpFileName);
214 else /* if (sumpParams.numRuns > 2) */
216 MrBayesPrint ("%s Summarizing parameters in %d files (%s.run1.p,\n", spacer, sumpParams.numRuns, sumpParams.sumpFileName);
217 MrBayesPrint ("%s %s.run2.p, etc)\n", spacer, sumpParams.sumpFileName);
219 MrBayesPrint ("%s Writing summary statistics to file %s.pstat\n", spacer, sumpParams.sumpFileName);
221 if (chainParams.relativeBurnin == YES)
222 MrBayesPrint ("%s Using relative burnin ('relburnin=yes'), discarding the first %.0f %% of samples\n",
223 spacer, chainParams.burninFraction*100.0, chainParams.burninFraction);
225 MrBayesPrint ("%s Using absolute burnin ('relburnin=no'), discarding the first %d samples\n",
226 spacer, chainParams.chainBurnIn, chainParams.chainBurnIn);
228 /* Initialize to silence warning. */
229 firstFileInfo.numRows = 0;
230 firstFileInfo.numColumns = 0;
232 /* examine input file(s) */
233 for (i=0; i<sumpParams.numRuns; i++)
235 if (sumpParams.numRuns == 1)
236 sprintf (temp, "%s.p", sumpParams.sumpFileName);
238 sprintf (temp, "%s.run%d.p", sumpParams.sumpFileName, i+1);
240 if (ExamineSumpFile (temp, &fileInfo, &headerNames, &nHeaders) == ERROR)
245 if (fileInfo.numRows == 0 || fileInfo.numColumns == 0)
247 MrBayesPrint ("%s The number of rows or columns in file %d is equal to zero\n", spacer, temp);
250 firstFileInfo = fileInfo;
254 if (firstFileInfo.numRows != fileInfo.numRows || firstFileInfo.numColumns != fileInfo.numColumns)
256 MrBayesPrint ("%s First file had %d rows and %d columns while file %s had %d rows and %d columns\n",
257 spacer, firstFileInfo.numRows, firstFileInfo.numColumns, temp, fileInfo.numRows, fileInfo.numColumns);
258 MrBayesPrint ("%s MrBayes expects the same number of rows and columns in all files\n", spacer);
264 numRows = fileInfo.numRows;
265 numColumns = fileInfo.numColumns;
266 numRuns = sumpParams.numRuns;
268 /* allocate space to hold parameter information */
269 if (AllocateParameterSamples (¶meterSamples, numRuns, numRows, numColumns) == ERROR)
273 for (i=0; i<sumpParams.numRuns; i++)
275 /* derive file name */
276 if (sumpParams.numRuns == 1)
277 sprintf (temp, "%s.p", sumpParams.sumpFileName);
279 sprintf (temp, "%s.run%d.p", sumpParams.sumpFileName, i+1);
282 if (ReadParamSamples (temp, &fileInfo, parameterSamples, i) == ERROR)
286 /* get length of longest header */
287 longestHeader = 9; /* 9 is the length of the word "parameter" (for printing table) */
288 for (i=0; i<nHeaders; i++)
290 len = (int) strlen(headerNames[i]);
291 if (len > longestHeader)
298 /* Print trace plots */
299 if (FindHeader("Gen", headerNames, nHeaders, &whichIsX) == ERROR)
301 MrBayesPrint ("%s Could not find the 'Gen' column\n", spacer);
304 if (FindHeader("LnL", headerNames, nHeaders, &whichIsY) == ERROR)
306 MrBayesPrint ("%s Could not find the 'LnL' column\n", spacer);
310 if (sumpParams.numRuns > 1)
312 if (sumpParams.allRuns == YES)
314 for (i=0; i<sumpParams.numRuns; i++)
316 MrBayesPrint ("\n%s Samples from run %d:\n", spacer, i+1);
317 if (PrintPlot (parameterSamples[whichIsX].values[i], parameterSamples[whichIsY].values[i], numRows) == ERROR)
323 if (PrintOverlayPlot (parameterSamples[whichIsX].values, parameterSamples[whichIsY].values, numRuns, 0, numRows) == ERROR)
329 if (PrintPlot (parameterSamples[whichIsX].values[0], parameterSamples[whichIsY].values[0], numRows) == ERROR)
333 /* calculate arithmetic and harmonic means of likelihoods */
335 /* open output file */
336 strncpy (temp, sumpParams.sumpOutfile, 90);
337 strcat (temp, ".lstat");
338 fpLstat = OpenNewMBPrintFile (temp);
342 /* print unique identifier to the output file */
343 if (strlen(stamp) > 1)
344 MrBayesPrintf (fpLstat, "[ID: %s]\n", stamp);
347 if (sumpParams.numRuns == 1)
348 MrBayesPrintf (fpLstat, "arithmetic_mean\tharmonic_mean\tvalues_discarded\n");
350 MrBayesPrintf (fpLstat, "run\tarithmetic_mean\tharmonic_mean\tvalues_discarded\n");
353 for (n=0; n<sumpParams.numRuns; n++)
356 if (HarmonicArithmeticMeanOnLogs (parameterSamples[whichIsY].values[n], numRows, &mean, &harm_mean) == ERROR)
361 if (sumpParams.numRuns == 1)
364 MrBayesPrint ("%s Estimated marginal likelihoods for run sampled in file \"%s.p\":\n", spacer, sumpParams.sumpFileName);
365 MrBayesPrint ("%s (Use the harmonic mean for Bayes factor comparisons of models)\n", spacer, sumpParams.sumpFileName);
366 MrBayesPrint ("%s (Values are saved to the file %s.lstat)\n\n", spacer, sumpParams.sumpOutfile);
367 MrBayesPrint ("%s Arithmetic mean Harmonic mean\n", spacer);
368 MrBayesPrint ("%s --------------------------------\n", spacer);
369 if (unreliable == NO)
370 MrBayesPrint ("%s %9.2lf %9.2lf\n", spacer, mean, harm_mean);
372 MrBayesPrint ("%s %9.2lf * %9.2lf *\n", spacer, mean, harm_mean);
375 MrBayesPrintf (fpLstat, "%s\t", MbPrintNum(mean));
376 MrBayesPrintf (fpLstat, "%s\t", MbPrintNum(harm_mean));
377 if (unreliable == YES)
378 MrBayesPrintf (fpLstat, "yes\n");
380 MrBayesPrintf (fpLstat, "no\n");
387 MrBayesPrint ("%s Estimated marginal likelihoods for runs sampled in files\n", spacer);
388 if (sumpParams.numRuns > 2)
389 MrBayesPrint ("%s \"%s.run1.p\", \"%s.run2.p\", etc:\n", spacer, sumpParams.sumpFileName, sumpParams.sumpFileName);
390 else /* if (sumpParams.numRuns == 2) */
391 MrBayesPrint ("%s \"%s.run1.p\" and \"%s.run2.p\":\n", spacer, sumpParams.sumpFileName, sumpParams.sumpFileName);
392 MrBayesPrint ("%s (Use the harmonic mean for Bayes factor comparisons of models)\n\n", spacer, sumpParams.sumpFileName);
393 MrBayesPrint ("%s (Values are saved to the file %s.lstat)\n\n", spacer, sumpParams.sumpOutfile);
394 MrBayesPrint ("%s Run Arithmetic mean Harmonic mean\n", spacer);
395 MrBayesPrint ("%s --------------------------------------\n", spacer);
397 if (unreliable == NO)
398 MrBayesPrint ("%s %3d %9.2lf %9.2lf\n", spacer, n+1, mean, harm_mean);
400 MrBayesPrint ("%s %3d %9.2lf * %9.2lf *\n", spacer, n+1, mean, harm_mean);
403 MrBayesPrintf (fpLstat, "%d\t", n+1);
404 MrBayesPrintf (fpLstat, "%s\t", MbPrintNum(mean));
405 MrBayesPrintf (fpLstat, "%s\t", MbPrintNum(harm_mean));
406 if (unreliable == YES)
407 MrBayesPrintf (fpLstat, "yes\n");
409 MrBayesPrintf (fpLstat, "no\n");
412 if (sumpParams.numRuns == 1)
414 MrBayesPrint ("%s --------------------------------\n", spacer);
418 if (HarmonicArithmeticMeanOnLogs (parameterSamples[whichIsY].values[0], sumpParams.numRuns*numRows, &mean, &harm_mean) == ERROR)
425 MrBayesPrint ("%s --------------------------------------\n", spacer);
426 if (unreliable == YES)
427 MrBayesPrint ("%s TOTAL %9.2lf * %9.2lf *\n", spacer, mean, harm_mean);
429 MrBayesPrint ("%s TOTAL %9.2lf %9.2lf\n", spacer, mean, harm_mean);
430 MrBayesPrint ("%s --------------------------------------\n", spacer);
432 /* print total to file */
433 MrBayesPrintf (fpLstat, "all\t");
434 MrBayesPrintf (fpLstat, "%s\t", MbPrintNum(mean));
435 MrBayesPrintf (fpLstat, "%s\t", MbPrintNum(harm_mean));
436 if (unreliable == YES)
437 MrBayesPrintf (fpLstat, "yes\n");
439 MrBayesPrintf (fpLstat, "no\n");
441 if (oneUnreliable == YES)
443 MrBayesPrint ("%s * These estimates may be unreliable because \n", spacer);
444 MrBayesPrint ("%s some extreme values were excluded\n\n", spacer);
450 SafeFclose(&fpLstat);
452 /* Calculate burnin */
453 burnin = fileInfo.firstParamLine - fileInfo.headerLine;
455 /* Print parameter information to screen and to file. */
456 if (sumpParams.numRuns > 1 && sumpParams.allRuns == YES)
458 for (i=0; i<sumpParams.numRuns; i++)
460 /* print table header */
462 MrBayesPrint ("%s Model parameter summaries for run sampled in file \"%s.run%d.p\":\n", spacer, sumpParams.sumpFileName, i+1);
463 MrBayesPrint ("%s (Based on %d samples out of a total of %d samples from this analysis)\n\n", spacer, numRows, numRows + burnin);
464 if (PrintParamStats (sumpParams.sumpOutfile, headerNames, nHeaders, parameterSamples, numRuns, numRows) == ERROR)
466 if (PrintModelStats (sumpParams.sumpOutfile, headerNames, nHeaders, parameterSamples, numRuns, numRows) == ERROR)
472 if (sumpParams.numRuns == 1)
473 MrBayesPrint ("%s Model parameter summaries for run sampled in file \"%s\":\n", spacer, sumpParams.sumpFileName);
474 else if (sumpParams.numRuns == 2)
476 MrBayesPrint ("%s Model parameter summaries over the runs sampled in files\n", spacer);
477 MrBayesPrint ("%s \"%s.run1.p\" and \"%s.run2.p\":\n", spacer, sumpParams.sumpFileName, sumpParams.sumpFileName);
481 MrBayesPrint ("%s Model parameter summaries over all %d runs sampled in files\n", spacer, sumpParams.numRuns);
482 MrBayesPrint ("%s \"%s.run1.p\", \"%s.run2.p\" etc:\n", spacer, sumpParams.sumpFileName, sumpParams.sumpFileName);
485 if (sumpParams.numRuns == 1)
487 MrBayesPrint ("%s Based on a total of %d samples out of a total of %d samples\n", spacer, numRows, numRows + burnin);
488 MrBayesPrint ("%s from this analysis.\n", spacer);
492 MrBayesPrint ("%s Summaries are based on a total of %d samples from %d runs.\n", spacer, sumpParams.numRuns*numRows, sumpParams.numRuns);
493 MrBayesPrint ("%s Each run produced %d samples of which %d samples were included.\n", spacer, numRows + burnin, numRows);
495 MrBayesPrint ("%s Parameter summaries saved to file \"%s.pstat\".\n", spacer, sumpParams.sumpOutfile);
497 if (PrintParamStats (sumpParams.sumpOutfile, headerNames, nHeaders, parameterSamples, numRuns, numRows) == ERROR)
499 if (PrintModelStats (sumpParams.sumpOutfile, headerNames, nHeaders, parameterSamples, numRuns, numRows) == ERROR)
503 FreeParameterSamples(parameterSamples);
504 for (i=0; i<nHeaders; i++)
505 free (headerNames[i]);
508 expecting = Expecting(COMMAND);
516 FreeParameterSamples (parameterSamples);
517 for (i=0; i<nHeaders; i++)
518 free (headerNames[i]);
522 SafeFclose (&fpLstat);
524 expecting = Expecting(COMMAND);
533 int i, nHeaders=0, numRows, numColumns, numRuns, whichIsX, whichIsY,
535 char **headerNames=NULL, temp[120];
536 SumpFileInfo fileInfo, firstFileInfo;
537 ParameterSample *parameterSamples=NULL;
538 int stepIndexSS,numSamplesInStepSS, stepBeginSS, stepBurnin;
539 MrBFlt *lnlp, *nextSteplnlp, *firstlnlp;
540 MrBFlt *marginalLnLSS=NULL,stepScalerSS,stepAcumulatorSS, stepLengthSS, tmpMfl;
541 int beginPrint, countPrint;
543 MrBFlt **plotArrayY=NULL,**plotArrayX=NULL;
548 # if defined (MPI_ENABLED)
553 chainParams.isSS=YES;
555 /* tell user we are ready to go */
556 if (sumssParams.numRuns == 1)
557 MrBayesPrint ("%s Summarizing parameters in file %s.p\n", spacer, sumpParams.sumpFileName);
558 else if (sumssParams.numRuns == 2)
559 MrBayesPrint ("%s Summarizing parameters in files %s.run1.p and %s.run2.p\n", spacer, sumpParams.sumpFileName, sumpParams.sumpFileName);
560 else /* if (sumssParams.numRuns > 2) */
562 MrBayesPrint ("%s Summarizing parameters in %d files (%s.run1.p,\n", spacer, sumssParams.numRuns, sumpParams.sumpFileName);
563 MrBayesPrint ("%s %s.run2.p, etc)\n", spacer, sumpParams.sumpFileName);
565 //MrBayesPrint ("%s Writing summary statistics to file %s.pstat\n", spacer, sumpParams.sumpFileName);
567 if (chainParams.relativeBurnin == YES)
568 MrBayesPrint ("%s Using relative burnin ('relburnin=yes'), discarding the first %.0f %% of samples within each step.\n",
569 spacer, chainParams.burninFraction*100.0, chainParams.burninFraction);
571 MrBayesPrint ("%s Using absolute burnin ('relburnin=no'), discarding the first %d samples within each step.\n",
572 spacer, chainParams.chainBurnIn, chainParams.chainBurnIn);
574 /* Initialize to silence warning. */
575 firstFileInfo.numRows = 0;
576 firstFileInfo.numColumns = 0;
578 /* examine input file(s) */
579 for (i=0; i<sumssParams.numRuns; i++)
581 if (sumssParams.numRuns == 1)
582 sprintf (temp, "%s.p", sumpParams.sumpFileName);
584 sprintf (temp, "%s.run%d.p", sumpParams.sumpFileName, i+1);
586 if (ExamineSumpFile (temp, &fileInfo, &headerNames, &nHeaders) == ERROR)
591 if (fileInfo.numRows == 0 || fileInfo.numColumns == 0)
593 MrBayesPrint ("%s The number of rows or columns in file %d is equal to zero\n", spacer, temp);
596 firstFileInfo = fileInfo;
600 if (firstFileInfo.numRows != fileInfo.numRows || firstFileInfo.numColumns != fileInfo.numColumns)
602 MrBayesPrint ("%s First file had %d rows and %d columns while file %s had %d rows and %d columns\n",
603 spacer, firstFileInfo.numRows, firstFileInfo.numColumns, temp, fileInfo.numRows, fileInfo.numColumns);
604 MrBayesPrint ("%s MrBayes expects the same number of rows and columns in all files\n", spacer);
610 numRows = fileInfo.numRows;
611 numColumns = fileInfo.numColumns;
612 numRuns = sumssParams.numRuns;
614 /* allocate space to hold parameter information */
615 if (AllocateParameterSamples (¶meterSamples, numRuns, numRows, numColumns) == ERROR)
619 for (i=0; i<sumssParams.numRuns; i++)
621 /* derive file name */
622 if (sumssParams.numRuns == 1)
623 sprintf (temp, "%s.p", sumpParams.sumpFileName);
625 sprintf (temp, "%s.run%d.p", sumpParams.sumpFileName, i+1);
628 if (ReadParamSamples (temp, &fileInfo, parameterSamples, i) == ERROR)
632 /* get length of longest header */
633 longestHeader = 9; /* 9 is the length of the word "parameter" (for printing table) */
634 for (i=0; i<nHeaders; i++)
636 len = (int) strlen(headerNames[i]);
637 if (len > longestHeader)
641 /* Print trace plots */
642 if (FindHeader("Gen", headerNames, nHeaders, &whichIsX) == ERROR)
644 MrBayesPrint ("%s Could not find the 'Gen' column\n", spacer);
647 if (FindHeader("LnL", headerNames, nHeaders, &whichIsY) == ERROR)
649 MrBayesPrint ("%s Could not find the 'LnL' column\n", spacer);
653 if (chainParams.burninSS > 0)
655 stepBeginSS = chainParams.burninSS + 1;
659 numSamplesInStepSS = (numRows-1)/(chainParams.numStepsSS-chainParams.burninSS);
660 stepBeginSS = numSamplesInStepSS + 1;
663 numSamplesInStepSS = (numRows - stepBeginSS)/chainParams.numStepsSS;
664 if ((numRows - stepBeginSS)%chainParams.numStepsSS!=0)
666 MrBayesPrint ("%s Error: Number of samples could not be evenly devided among steps (%d samples among %d steps). \n", spacer,(numRows - stepBeginSS),chainParams.numStepsSS);
670 if (chainParams.relativeBurnin == YES)
672 stepBurnin = (int)(numSamplesInStepSS*chainParams.burninFraction);
676 stepBurnin = chainParams.chainBurnIn;
677 if (stepBurnin >= numSamplesInStepSS)
679 MrBayesPrint ("%s Error: Burnin in each step(%d) is longer then the step itself(%d). \n", spacer,stepBurnin, numSamplesInStepSS);
684 marginalLnLSS = (MrBFlt *) SafeCalloc (sumssParams.numRuns, sizeof(MrBFlt));
685 /*Preparing and printing joined plot.*/
686 plotArrayY = (MrBFlt **) SafeCalloc (sumssParams.numRuns+1, sizeof(MrBFlt*));
687 for (i=0; i<sumssParams.numRuns+1; i++)
688 plotArrayY[i] = (MrBFlt *) SafeCalloc (numSamplesInStepSS, sizeof(MrBFlt));
690 plotArrayX = (MrBFlt **) SafeCalloc (sumssParams.numRuns, sizeof(MrBFlt*));
691 for (i=0; i<sumssParams.numRuns; i++)
693 plotArrayX[i] = (MrBFlt *) SafeCalloc (numSamplesInStepSS, sizeof(MrBFlt));
694 for (j=0; j<numSamplesInStepSS; j++)
695 plotArrayX[i][j]=j+1;
698 MrBayesPrint ("%s In total %d sampls are red from .p files.\n", spacer, numRows);
700 MrBayesPrint ("%s Marginal likelihood (in natural log units) is estimated using stepping-stone sampling\n", spacer);
701 MrBayesPrint ("%s based on %d steps with %d samples within each step. \n", spacer, chainParams.numStepsSS, numSamplesInStepSS);
702 MrBayesPrint ("%s First %d samples (including generation 0) are discarded as initial burn-in.\n", spacer, stepBeginSS);
703 if (chainParams.startFromPriorSS==YES)
704 MrBayesPrint ("%s Sampling is assumed have being done from prior to posterior.\n", spacer);
707 MrBayesPrint ("%s Sampling is assumed have being done from posterior to prior.\n", spacer);
712 MrBayesPrint ("\n\n%s Step contribution table.\n\n",spacer);
713 MrBayesPrint (" Columns in the table: \n");
714 MrBayesPrint (" Step -- Index of the step \n");
715 MrBayesPrint (" runX -- Contribution to the marginal log likelihood of run X, i.e. marginal \n");
716 MrBayesPrint (" log likelihood for run X is the sum across all steps in column runX.\n\n");
718 if (firstPass == YES && chainParams.relativeBurnin == YES)
719 MrBayesPrint ("%s The table entrances are based on samples excluding burn-in %d samples (%d%%) \n", spacer, stepBurnin,(int)(100*chainParams.burninFraction));
721 MrBayesPrint ("%s The table entrances are based on samples excluding burn-in %d samples \n", spacer, stepBurnin);
722 MrBayesPrint ("%s discarded at the begining of each step. \n\n", spacer);
724 //MrBayesPrint ("%s Run Marginal likelihood (ln)\n",spacer);
725 //MrBayesPrint ("%s ------------------------------\n",spacer);
726 MrBayesPrint (" Step");
727 for (j=0; j<sumssParams.numRuns ; j++)
731 MrBayesPrint (" run%d", j+1);
734 for (i=0; i<sumssParams.numRuns; i++)
736 marginalLnLSS[i] = 0.0;
738 for (stepIndexSS = chainParams.numStepsSS-1; stepIndexSS>=0; stepIndexSS--)
740 if (chainParams.startFromPriorSS==YES)
742 stepLengthSS = BetaQuantile(chainParams.alphaSS, 1.0, (MrBFlt)(chainParams.numStepsSS-stepIndexSS)/(MrBFlt)chainParams.numStepsSS) - BetaQuantile(chainParams.alphaSS, 1.0, (MrBFlt)(chainParams.numStepsSS-1-stepIndexSS)/(MrBFlt)chainParams.numStepsSS);
746 stepLengthSS = BetaQuantile(chainParams.alphaSS, 1.0, (MrBFlt)(stepIndexSS+1)/(MrBFlt)chainParams.numStepsSS) - BetaQuantile(chainParams.alphaSS, 1.0, (MrBFlt)stepIndexSS/(MrBFlt)chainParams.numStepsSS);
748 MrBayesPrint (" %3d ", chainParams.numStepsSS-stepIndexSS);
749 for (i=0; i<sumssParams.numRuns; i++)
751 lnlp = parameterSamples[whichIsY].values[i] + stepBeginSS + (chainParams.numStepsSS-stepIndexSS-1)*numSamplesInStepSS;
752 nextSteplnlp = lnlp+numSamplesInStepSS;
754 stepAcumulatorSS = 0.0;
755 stepScalerSS = *lnlp*stepLengthSS;
756 while (lnlp < nextSteplnlp)
758 if (*lnlp*stepLengthSS > stepScalerSS + 200.0)
761 stepAcumulatorSS /= exp(*lnlp*stepLengthSS - 100.0 - stepScalerSS);
762 stepScalerSS= *lnlp*stepLengthSS - 100.0;
764 stepAcumulatorSS += exp(*lnlp*stepLengthSS - stepScalerSS);
767 tmpMfl = (log(stepAcumulatorSS/(numSamplesInStepSS-stepBurnin)) + stepScalerSS);
768 MrBayesPrint (" %10.3lf", tmpMfl);
769 marginalLnLSS[i] += tmpMfl;
772 //MrBayesPrint ("%s %3d %9.2f \n", spacer, i+1, marginalLnLSS);
775 for (j=0; j<sumssParams.numRuns ; j++)
779 MrBayesPrint ("----------");
782 MrBayesPrint (" Sum: ");
783 for (j=0; j<sumssParams.numRuns ; j++)
784 MrBayesPrint (" %10.3lf", marginalLnLSS[j]);
789 if (sumssParams.numRuns > 1)
791 MrBayesPrint ("%s Below are rough plots of the generations (x-axis) during burn in \n", spacer);
792 MrBayesPrint ("%s phase versus the log probability of observing the data (y-axis). \n", spacer);
793 MrBayesPrint ("%s You can use these graphs to determine if the burn in for your SS \n", spacer);
794 MrBayesPrint ("%s analysis was sufficiant. The log probability suppose to plateau \n", spacer);
795 MrBayesPrint ("%s indicating that you may be at stationarity by the time you finish \n", spacer);
796 MrBayesPrint ("%s burn in phase. This burn in, unlike burn in within each step, is \n", spacer);
797 MrBayesPrint ("%s fixed and can not be changed. \n", spacer);
801 MrBayesPrint ("%s Below is a rough plot of the generations (x-axis) during burn in \n", spacer);
802 MrBayesPrint ("%s phase versus the log probability of observing the data (y-axis). \n", spacer);
803 MrBayesPrint ("%s You can use these graph to determine if the burn in for your SS \n", spacer);
804 MrBayesPrint ("%s analysis was sufficiant. The log probability suppose to plateau \n", spacer);
805 MrBayesPrint ("%s indicating that you may be at stationarity by the time you finish \n", spacer);
806 MrBayesPrint ("%s burn in phase. This burn in, unlike burn in within each step, is \n", spacer);
807 MrBayesPrint ("%s fixed and can not be changed. \n", spacer);
812 goto sumssExitOptions;
816 MrBayesPrint ("\n\n%s Step plot(s).\n",spacer);
820 if (sumssParams.stepToPlot == 0)
822 beginPrint=(int)(sumssParams.discardFraction*stepBeginSS);
823 countPrint=stepBeginSS-beginPrint;
824 MrBayesPrint ("%s Ploting step 0, i.e initial burn-in phase consisting of %d samples.\n", spacer,stepBeginSS);
825 MrBayesPrint ("%s According to 'Discardfrac=%.2f', first %d samples are not ploted.\n", spacer,sumssParams.discardFraction,beginPrint);
829 if (sumssParams.stepToPlot > chainParams.numStepsSS)
831 MrBayesPrint ("%s Chosen index of step to print %d is out of range of step indices[0,...,%d].\n", spacer,sumssParams.stepToPlot,chainParams.numStepsSS);
834 beginPrint=stepBeginSS+(sumssParams.stepToPlot-1)*numSamplesInStepSS + (int)(sumssParams.discardFraction*numSamplesInStepSS);
835 countPrint=numSamplesInStepSS-(int)(sumssParams.discardFraction*numSamplesInStepSS);
836 MrBayesPrint ("%s Ploting step %d consisting of %d samples.\n", spacer,sumssParams.stepToPlot,numSamplesInStepSS);
837 MrBayesPrint ("%s According to 'Discardfrac=%.2f', first %d samples are not ploted.\n", spacer,sumssParams.discardFraction,(int)(sumssParams.discardFraction*numSamplesInStepSS));
840 if (sumssParams.numRuns > 1)
842 if (sumpParams.allRuns == YES)
844 for (i=0; i<sumssParams.numRuns; i++)
846 MrBayesPrint ("\n%s Samples from run %d:\n", spacer, i+1);
847 if (PrintPlot (parameterSamples[whichIsX].values[i]+beginPrint, parameterSamples[whichIsY].values[i]+beginPrint, countPrint) == ERROR)
853 if (PrintOverlayPlot (parameterSamples[whichIsX].values, parameterSamples[whichIsY].values, numRuns, beginPrint, countPrint) == ERROR)
859 if (PrintPlot (parameterSamples[whichIsX].values[0]+beginPrint, parameterSamples[whichIsY].values[0]+beginPrint, countPrint) == ERROR)
863 if (sumssParams.askForMorePlots == NO || firstPass == YES)
866 MrBayesPrint (" You can choose to print new step plots for different steps or discard fractions.\n");
867 MrBayesPrint (" Allowed range of 'Steptoplot' are from 0 to %d.\n", chainParams.numStepsSS);
868 MrBayesPrint (" If the next entered value is negative, 'sumss' will stop printing step plots.\n");
869 MrBayesPrint (" If the next entered value is positive, but out of range, you will be offered\n");
870 MrBayesPrint (" to change parameter 'Discardfrac' of 'sumss'.\n");
871 MrBayesPrint (" Enter new step number 'Steptoplot':");
875 if (j > chainParams.numStepsSS)
879 MrBayesPrint (" Enter new value for 'Discardfrac', should be in range 0.0 to 1.0:");
880 k = scanf("%f",&tmpf);
881 sumssParams.discardFraction = (MrBFlt)tmpf;
883 while (sumssParams.discardFraction < 0.0 || sumssParams.discardFraction > 1.0);
886 sumssParams.stepToPlot=j;
890 goto sumssExitOptions;
894 MrBayesPrint ("\n\n%s Joined plot(s).\n",spacer);
898 MrBayesPrint ("%s Joined plot of %d samples of all steps together. 'smoothing' is set to:%d\n", spacer,numSamplesInStepSS,sumssParams.smoothing);
899 MrBayesPrint ("%s According to step burn-in, first %d samples are not ploted.\n", spacer,stepBurnin);
901 for (i=0; i<sumssParams.numRuns; i++)
903 for (j=stepBurnin;j<numSamplesInStepSS;j++)
904 plotArrayY[sumssParams.numRuns][j]=0.0;
905 lnlp= parameterSamples[whichIsY].values[i] + stepBeginSS;
907 for (stepIndexSS = chainParams.numStepsSS-1; stepIndexSS>0; stepIndexSS--)
909 firstlnlp=plotArrayY[sumssParams.numRuns] + stepBurnin;
911 nextSteplnlp += numSamplesInStepSS;
912 while (lnlp < nextSteplnlp)
919 for (j=stepBurnin;j<numSamplesInStepSS;j++)
923 for (k=j-sumssParams.smoothing;k<=j+sumssParams.smoothing;k++)
925 if (k>=stepBurnin && k<numSamplesInStepSS)
927 sum += plotArrayY[sumssParams.numRuns][k];
931 plotArrayY[i][j] = sum/count;
933 if (max < plotArrayY[i][j])
934 max=plotArrayY[i][j];
937 /* for (j=stepBurnin;j<numSamplesInStepSS;j++)
939 plotArrayY[i][j] /= max;
943 beginPrint=stepBurnin;
944 countPrint=numSamplesInStepSS-stepBurnin;
946 if (sumssParams.numRuns > 1)
948 if (sumpParams.allRuns == YES)
950 for (i=0; i<sumssParams.numRuns; i++)
952 MrBayesPrint ("\n%s Samples from run %d:\n", spacer, i+1);
953 if (PrintPlot (plotArrayX[i]+beginPrint, plotArrayY[i]+beginPrint, countPrint) == ERROR)
959 if (PrintOverlayPlot (plotArrayX, plotArrayY, numRuns, beginPrint, countPrint) == ERROR)
965 if (PrintPlot (plotArrayX[0]+beginPrint, plotArrayY[0]+beginPrint, countPrint) == ERROR)
969 if (sumssParams.askForMorePlots == NO || firstPass == YES)
972 MrBayesPrint (" You can choose to print new joined plots with different step burn-in or smoothing.\n");
973 MrBayesPrint (" Allowed range of step burn-in values are from 0 to %d.\n", numSamplesInStepSS-1);
974 MrBayesPrint (" If the next entered value is negative, 'sumss' will stop printing joined plots.\n");
975 MrBayesPrint (" If the next entered value is positive, but out of range, you will be offered\n");
976 MrBayesPrint (" to change 'Smoothing'.\n");
977 MrBayesPrint (" Enter new step burn-in:");
981 if (j >= numSamplesInStepSS)
983 MrBayesPrint (" Enter new value for 'Smoothing':");
985 sumssParams.smoothing = abs(j);
995 if (sumssParams.askForMorePlots == YES)
998 MrBayesPrint (" Sumss is interactive, because of parameter 'Askmore=YES' setting. \n");
999 MrBayesPrint (" What would you like to do next?\n");
1000 MrBayesPrint (" 1) Print updated table according to new step burn-in.\n");
1001 MrBayesPrint (" 2) Print Step plot(s).\n");
1002 MrBayesPrint (" 3) Print Joined plot(s).\n");
1003 MrBayesPrint (" 4) Exit 'sumss'.\n");
1004 MrBayesPrint (" Enter a number that corresponds to one of the options:");
1008 }while (j<1 || j>4);
1012 MrBayesPrint (" Allowed range of step burn-in values are from 0 to %d\n", numSamplesInStepSS-1);
1013 MrBayesPrint (" Current step burn-in value is:%d\n", stepBurnin);
1014 MrBayesPrint (" Enter new step burn-in:");
1017 k = scanf("%d",&stepBurnin);
1019 while (stepBurnin < 0 || stepBurnin > numSamplesInStepSS-1);
1020 MrBayesPrint ("\n");
1028 goto sumssJoinedPlot;
1032 FreeParameterSamples(parameterSamples);
1033 for (i=0; i<nHeaders; i++)
1034 free (headerNames[i]);
1037 expecting = Expecting(COMMAND);
1038 strcpy (spacer, "");
1039 chainParams.isSS=NO;
1040 for (i=0; i<sumssParams.numRuns+1; i++)
1041 free(plotArrayY[i]);
1043 for (i=0; i<sumssParams.numRuns; i++)
1044 free(plotArrayX[i]);
1046 free(marginalLnLSS);
1053 FreeParameterSamples (parameterSamples);
1054 if (headerNames!=NULL)
1055 for (i=0; i<nHeaders; i++)
1056 free (headerNames[i]);
1059 expecting = Expecting(COMMAND);
1060 strcpy (spacer, "");
1061 chainParams.isSS=NO;
1062 if (plotArrayY!=NULL)
1063 for (i=0; i<sumssParams.numRuns+1; i++)
1064 free(plotArrayY[i]);
1066 if (plotArrayX!=NULL)
1067 for (i=0; i<sumssParams.numRuns; i++)
1068 free(plotArrayX[i]);
1070 free(marginalLnLSS);
1076 int DoSumpParm (char *parmName, char *tkn)
1082 if (expecting == Expecting(PARAMETER))
1084 expecting = Expecting(EQUALSIGN);
1088 if (!strcmp(parmName, "Xxxxxxxxxx"))
1090 expecting = Expecting(PARAMETER);
1091 expecting |= Expecting(SEMICOLON);
1093 /* set Filename (sumpParams.sumpFileName) ***************************************************/
1094 else if (!strcmp(parmName, "Filename"))
1096 if (expecting == Expecting(EQUALSIGN))
1098 expecting = Expecting(ALPHA);
1101 else if (expecting == Expecting(ALPHA))
1103 if (strlen(tkn)>99 && (strchr(tkn,' ')-tkn) > 99)
1105 MrBayesPrint ("%s Maximum allowed length of file name is 99 characters. The given name:\n", spacer);
1106 MrBayesPrint ("%s '%s'\n", spacer,tkn);
1109 sscanf (tkn, "%s", tempStr);
1110 strcpy (sumpParams.sumpFileName, tempStr);
1111 strcpy (sumpParams.sumpOutfile, tempStr);
1112 MrBayesPrint ("%s Setting sump filename and output file name to %s\n", spacer, sumpParams.sumpFileName);
1113 expecting = Expecting(PARAMETER) | Expecting(SEMICOLON);
1118 /* set Outputname (sumpParams.sumpOutfile) *******************************************************/
1119 else if (!strcmp(parmName, "Outputname"))
1121 if (expecting == Expecting(EQUALSIGN))
1123 expecting = Expecting(ALPHA);
1126 else if (expecting == Expecting(ALPHA))
1128 if (strlen(tkn)>99 && (strchr(tkn,' ')-tkn) > 99)
1130 MrBayesPrint ("%s Maximum allowed length of file name is 99 characters. The given name:\n", spacer);
1131 MrBayesPrint ("%s '%s'\n", spacer,tkn);
1134 sscanf (tkn, "%s", tempStr);
1135 strcpy (sumpParams.sumpOutfile, tempStr);
1136 MrBayesPrint ("%s Setting sump output file name to \"%s\"\n", spacer, sumpParams.sumpOutfile);
1137 expecting = Expecting(PARAMETER) | Expecting(SEMICOLON);
1142 /* set Relburnin (chainParams.relativeBurnin) ********************************************************/
1143 else if (!strcmp(parmName, "Relburnin"))
1145 if (expecting == Expecting(EQUALSIGN))
1146 expecting = Expecting(ALPHA);
1147 else if (expecting == Expecting(ALPHA))
1149 if (IsArgValid(tkn, tempStr) == NO_ERROR)
1151 if (!strcmp(tempStr, "Yes"))
1152 chainParams.relativeBurnin = YES;
1154 chainParams.relativeBurnin = NO;
1158 MrBayesPrint ("%s Invalid argument for Relburnin\n", spacer);
1161 if (chainParams.relativeBurnin == YES)
1162 MrBayesPrint ("%s Using relative burnin (a fraction of samples discarded).\n", spacer);
1164 MrBayesPrint ("%s Using absolute burnin (a fixed number of samples discarded).\n", spacer);
1165 expecting = Expecting(PARAMETER) | Expecting(SEMICOLON);
1172 /* set Burnin (chainParams.chainBurnIn) ***********************************************************/
1173 else if (!strcmp(parmName, "Burnin"))
1175 if (expecting == Expecting(EQUALSIGN))
1176 expecting = Expecting(NUMBER);
1177 else if (expecting == Expecting(NUMBER))
1179 sscanf (tkn, "%d", &tempI);
1180 chainParams.chainBurnIn = tempI;
1181 MrBayesPrint ("%s Setting burn-in to %d\n", spacer, chainParams.chainBurnIn);
1182 expecting = Expecting(PARAMETER) | Expecting(SEMICOLON);
1189 /* set Burninfrac (chainParams.burninFraction) ************************************************************/
1190 else if (!strcmp(parmName, "Burninfrac"))
1192 if (expecting == Expecting(EQUALSIGN))
1193 expecting = Expecting(NUMBER);
1194 else if (expecting == Expecting(NUMBER))
1196 sscanf (tkn, "%lf", &tempD);
1199 MrBayesPrint ("%s Burnin fraction too low (< 0.01)\n", spacer);
1204 MrBayesPrint ("%s Burnin fraction too high (> 0.50)\n", spacer);
1207 chainParams.burninFraction = tempD;
1208 MrBayesPrint ("%s Setting burnin fraction to %.2f\n", spacer, chainParams.burninFraction);
1209 expecting = Expecting(PARAMETER) | Expecting(SEMICOLON);
1216 /* set Minprob (sumpParams.minProb) ************************************************************/
1217 else if (!strcmp(parmName, "Minprob"))
1219 if (expecting == Expecting(EQUALSIGN))
1220 expecting = Expecting(NUMBER);
1221 else if (expecting == Expecting(NUMBER))
1223 sscanf (tkn, "%lf", &tempD);
1226 MrBayesPrint ("%s Minprob too high (it should be smaller than 0.50)\n", spacer);
1229 sumpParams.minProb = tempD;
1230 MrBayesPrint ("%s Setting minprob to %1.3f\n", spacer, sumpParams.minProb);
1231 expecting = Expecting(PARAMETER) | Expecting(SEMICOLON);
1238 /* set Nruns (sumpParams.numRuns) *******************************************************/
1239 else if (!strcmp(parmName, "Nruns"))
1241 if (expecting == Expecting(EQUALSIGN))
1242 expecting = Expecting(NUMBER);
1243 else if (expecting == Expecting(NUMBER))
1245 sscanf (tkn, "%d", &tempI);
1248 MrBayesPrint ("%s Nruns must be at least 1\n", spacer);
1253 sumpParams.numRuns = tempI;
1254 MrBayesPrint ("%s Setting sump nruns to %d\n", spacer, sumpParams.numRuns);
1255 expecting = Expecting(PARAMETER) | Expecting(SEMICOLON);
1261 /* set Hpd (sumpParams.HPD) ********************************************************/
1262 else if (!strcmp(parmName, "Hpd"))
1264 if (expecting == Expecting(EQUALSIGN))
1265 expecting = Expecting(ALPHA);
1266 else if (expecting == Expecting(ALPHA))
1268 if (IsArgValid(tkn, tempStr) == NO_ERROR)
1270 if (!strcmp(tempStr, "Yes"))
1271 sumpParams.HPD = YES;
1273 sumpParams.HPD = NO;
1277 MrBayesPrint ("%s Invalid argument for Hpd\n", spacer);
1280 if (sumpParams.HPD == YES)
1281 MrBayesPrint ("%s Reporting 95 %% region of Highest Posterior Density (HPD).\n", spacer);
1283 MrBayesPrint ("%s Reporting median interval containing 95 %% of values.\n", spacer);
1284 expecting = Expecting(PARAMETER) | Expecting(SEMICOLON);
1291 /* set Allruns (sumpParams.allRuns) ********************************************************/
1292 else if (!strcmp(parmName, "Allruns"))
1294 if (expecting == Expecting(EQUALSIGN))
1295 expecting = Expecting(ALPHA);
1296 else if (expecting == Expecting(ALPHA))
1298 if (IsArgValid(tkn, tempStr) == NO_ERROR)
1300 if (!strcmp(tempStr, "Yes"))
1301 sumpParams.allRuns = YES;
1303 sumpParams.allRuns = NO;
1307 MrBayesPrint ("%s Invalid argument for allruns (valid arguments are 'yes' and 'no')\n", spacer);
1310 if (sumpParams.allRuns == YES)
1311 MrBayesPrint ("%s Setting sump to print information for each run\n", spacer);
1313 MrBayesPrint ("%s Setting sump to print only summary information for all runs\n", spacer);
1314 expecting = Expecting(PARAMETER) | Expecting(SEMICOLON);
1327 int DoSumSsParm (char *parmName, char *tkn)
1333 if (expecting == Expecting(PARAMETER))
1335 expecting = Expecting(EQUALSIGN);
1339 if (!strcmp(parmName, "Xxxxxxxxxx"))
1341 expecting = Expecting(PARAMETER);
1342 expecting |= Expecting(SEMICOLON);
1344 /* set Filename (sumpParams.sumpFileName) ***************************************************/
1345 else if (!strcmp(parmName, "Filename"))
1347 if (expecting == Expecting(EQUALSIGN))
1349 expecting = Expecting(ALPHA);
1352 else if (expecting == Expecting(ALPHA))
1354 if (strlen(tkn)>99 && (strchr(tkn,' ')-tkn) > 99)
1356 MrBayesPrint ("%s Maximum allowed length of file name is 99 characters. The given name:\n", spacer);
1357 MrBayesPrint ("%s '%s'\n", spacer,tkn);
1360 sscanf (tkn, "%s", tempStr);
1361 strcpy (sumpParams.sumpFileName, tempStr);
1362 strcpy (sumpParams.sumpOutfile, tempStr);
1363 MrBayesPrint ("%s Setting sump filename and output file name to %s\n", spacer, sumpParams.sumpFileName);
1364 expecting = Expecting(PARAMETER) | Expecting(SEMICOLON);
1369 /* set Outputname (sumpParams.sumpOutfile) *******************************************************/
1370 /*else if (!strcmp(parmName, "Outputname"))
1372 if (expecting == Expecting(EQUALSIGN))
1374 expecting = Expecting(ALPHA);
1377 else if (expecting == Expecting(ALPHA))
1379 if (strlen(tkn)>99 && (strchr(tkn,' ')-tkn) > 99)
1381 MrBayesPrint ("%s Maximum allowed length of file name is 99 characters. The given name:\n", spacer);
1382 MrBayesPrint ("%s '%s'\n", spacer,tkn);
1385 sscanf (tkn, "%s", tempStr);
1386 strcpy (sumpParams.sumpOutfile, tempStr);
1387 MrBayesPrint ("%s Setting sump output file name to \"%s\"\n", spacer, sumpParams.sumpOutfile);
1388 expecting = Expecting(PARAMETER) | Expecting(SEMICOLON);
1393 /* set Relburnin (chainParams.relativeBurnin) ********************************************************/
1394 else if (!strcmp(parmName, "Relburnin"))
1396 if (expecting == Expecting(EQUALSIGN))
1397 expecting = Expecting(ALPHA);
1398 else if (expecting == Expecting(ALPHA))
1400 if (IsArgValid(tkn, tempStr) == NO_ERROR)
1402 if (!strcmp(tempStr, "Yes"))
1403 chainParams.relativeBurnin = YES;
1405 chainParams.relativeBurnin = NO;
1409 MrBayesPrint ("%s Invalid argument for Relburnin\n", spacer);
1412 if (chainParams.relativeBurnin == YES)
1413 MrBayesPrint ("%s Using relative burnin (a fraction of samples discarded).\n", spacer);
1415 MrBayesPrint ("%s Using absolute burnin (a fixed number of samples discarded).\n", spacer);
1416 expecting = Expecting(PARAMETER) | Expecting(SEMICOLON);
1423 /* set Burnin (chainParams.chainBurnIn) ***********************************************************/
1424 else if (!strcmp(parmName, "Burnin"))
1426 if (expecting == Expecting(EQUALSIGN))
1427 expecting = Expecting(NUMBER);
1428 else if (expecting == Expecting(NUMBER))
1430 sscanf (tkn, "%d", &tempI);
1431 chainParams.chainBurnIn = tempI;
1432 MrBayesPrint ("%s Setting burn-in to %d\n", spacer, chainParams.chainBurnIn);
1433 expecting = Expecting(PARAMETER) | Expecting(SEMICOLON);
1440 /* set Burninfrac (chainParams.burninFraction) ************************************************************/
1441 else if (!strcmp(parmName, "Burninfrac"))
1443 if (expecting == Expecting(EQUALSIGN))
1444 expecting = Expecting(NUMBER);
1445 else if (expecting == Expecting(NUMBER))
1447 sscanf (tkn, "%lf", &tempD);
1450 MrBayesPrint ("%s Burnin fraction too low (< 0.01)\n", spacer);
1455 MrBayesPrint ("%s Burnin fraction too high (> 0.50)\n", spacer);
1458 chainParams.burninFraction = tempD;
1459 MrBayesPrint ("%s Setting burnin fraction to %.2f\n", spacer, chainParams.burninFraction);
1460 expecting = Expecting(PARAMETER) | Expecting(SEMICOLON);
1467 /* set Nruns (sumssParams.numRuns) *******************************************************/
1468 else if (!strcmp(parmName, "Nruns"))
1470 if (expecting == Expecting(EQUALSIGN))
1471 expecting = Expecting(NUMBER);
1472 else if (expecting == Expecting(NUMBER))
1474 sscanf (tkn, "%d", &tempI);
1477 MrBayesPrint ("%s Nruns must be at least 1\n", spacer);
1482 sumssParams.numRuns = tempI;
1483 MrBayesPrint ("%s Setting sumss nruns to %d\n", spacer, sumssParams.numRuns);
1484 expecting = Expecting(PARAMETER) | Expecting(SEMICOLON);
1490 /* set Allruns (sumssParams.allRuns) ********************************************************/
1491 else if (!strcmp(parmName, "Allruns"))
1493 if (expecting == Expecting(EQUALSIGN))
1494 expecting = Expecting(ALPHA);
1495 else if (expecting == Expecting(ALPHA))
1497 if (IsArgValid(tkn, tempStr) == NO_ERROR)
1499 if (!strcmp(tempStr, "Yes"))
1500 sumssParams.allRuns = YES;
1502 sumssParams.allRuns = NO;
1506 MrBayesPrint ("%s Invalid argument for allruns (valid arguments are 'yes' and 'no')\n", spacer);
1509 if (sumssParams.allRuns == YES)
1510 MrBayesPrint ("%s Setting sump to print information for each run\n", spacer);
1512 MrBayesPrint ("%s Setting sump to print only summary information for all runs\n", spacer);
1513 expecting = Expecting(PARAMETER) | Expecting(SEMICOLON);
1518 /* set Steptoplot (sumssParams.stepToPlot) *******************************************************/
1519 else if (!strcmp(parmName, "Steptoplot"))
1521 if (expecting == Expecting(EQUALSIGN))
1522 expecting = Expecting(NUMBER);
1523 else if (expecting == Expecting(NUMBER))
1525 sscanf (tkn, "%d", &tempI);
1528 MrBayesPrint ("%s Steptoplot must be at least 0\n", spacer);
1533 sumssParams.stepToPlot = tempI;
1534 MrBayesPrint ("%s Setting sumss steptoplot to %d\n", spacer, sumssParams.stepToPlot);
1535 expecting = Expecting(PARAMETER) | Expecting(SEMICOLON);
1541 /* set Smoothing (sumssParams.smoothing) *******************************************************/
1542 else if (!strcmp(parmName, "Smoothing"))
1544 if (expecting == Expecting(EQUALSIGN))
1545 expecting = Expecting(NUMBER);
1546 else if (expecting == Expecting(NUMBER))
1548 sscanf (tkn, "%d", &tempI);
1551 MrBayesPrint ("%s Smoothing must be at least 0\n", spacer);
1556 sumssParams.smoothing = tempI;
1557 MrBayesPrint ("%s Setting sumss smoothing to %d\n", spacer, sumssParams.smoothing);
1558 expecting = Expecting(PARAMETER) | Expecting(SEMICOLON);
1564 /* set Allruns (sumssParams.askForMorePlots) ********************************************************/
1565 else if (!strcmp(parmName, "Askmore"))
1567 if (expecting == Expecting(EQUALSIGN))
1568 expecting = Expecting(ALPHA);
1569 else if (expecting == Expecting(ALPHA))
1571 if (IsArgValid(tkn, tempStr) == NO_ERROR)
1573 if (!strcmp(tempStr, "Yes"))
1574 sumssParams.askForMorePlots = YES;
1576 sumssParams.askForMorePlots = NO;
1580 MrBayesPrint ("%s Invalid argument for askmore (valid arguments are 'yes' and 'no')\n", spacer);
1583 if (sumssParams.askForMorePlots == YES)
1584 MrBayesPrint ("%s Setting sumss to be interactiva by asking for more plots of burn-in or individual steps.\n", spacer);
1586 MrBayesPrint ("%s Setting sumss not to be interactive. It will not ask to print more plots.\n", spacer);
1587 expecting = Expecting(PARAMETER) | Expecting(SEMICOLON);
1592 /* set Discardfrac (sumssParams.discardFraction) ************************************************************/
1593 else if (!strcmp(parmName, "Discardfrac"))
1595 if (expecting == Expecting(EQUALSIGN))
1596 expecting = Expecting(NUMBER);
1597 else if (expecting == Expecting(NUMBER))
1599 sscanf (tkn, "%lf", &tempD);
1602 MrBayesPrint ("%s Discard fraction too low (< 0.00)\n", spacer);
1607 MrBayesPrint ("%s Discard fraction too high (> 1.00)\n", spacer);
1610 sumssParams.discardFraction = tempD;
1611 MrBayesPrint ("%s Setting discard fraction to %.2f\n", spacer, sumssParams.discardFraction);
1612 expecting = Expecting(PARAMETER) | Expecting(SEMICOLON);
1627 /* ExamineSumpFile: Collect info on the parameter samples in the file */
1628 int ExamineSumpFile (char *fileName, SumpFileInfo *fileInfo, char ***headerNames, int *nHeaders)
1630 char *sumpTokenP, sumpToken[CMD_STRING_LENGTH], *s=NULL, *headerLine, *t;
1631 int i, lineTerm, inSumpComment, lineNum, lastNonDigitLine, numParamLines, allDigitLine,
1632 lastTokenWasDash, nNumbersOnThisLine, tokenType, burnin, nLines, firstNumCols;
1636 /* open binary file */
1637 if ((fp = OpenBinaryFileR(fileName)) == NULL)
1639 /* test for some simple errors */
1640 if (strlen(sumpParams.sumpFileName) > 2)
1642 s = sumpParams.sumpFileName + (int) strlen(sumpParams.sumpFileName) - 2;
1643 if (strcmp(s, ".p") == 0)
1645 MrBayesPrint ("%s It appears that you need to remove '.p' from the 'Filename' parameter\n", spacer);
1646 MrBayesPrint ("%s Also make sure that 'Nruns' is set correctly\n", spacer);
1650 MrBayesPrint ("%s Make sure that 'Nruns' is set correctly\n", spacer);
1654 /* find out what type of line termination is used */
1655 lineTerm = LineTermType (fp);
1656 if (lineTerm != LINETERM_MAC && lineTerm != LINETERM_DOS && lineTerm != LINETERM_UNIX)
1658 MrBayesPrint ("%s Unknown line termination\n", spacer);
1662 /* find length of longest line */
1663 fileInfo->longestLineLength = LongestLine (fp);
1664 fileInfo->longestLineLength += 10; /* better safe than sorry; if you fgets with raw longestLineLength, you run into problems */
1666 /* allocate string long enough to hold a line */
1667 s = (char *)SafeMalloc(2 * ((size_t)(fileInfo->longestLineLength) +10) * sizeof(char));
1670 MrBayesPrint ("%s Problem allocating string for reading sump file\n", spacer);
1673 headerLine = s + fileInfo->longestLineLength + 10;
1675 /* close binary file */
1678 /* open text file */
1679 if ((fp = OpenTextFileR(fileName)) == NULL)
1682 /* Check file for appropriate blocks. We want to find the last block
1683 in the file and start from there. */
1685 lineNum = lastNonDigitLine = numParamLines = 0;
1686 while (fgets (s, fileInfo->longestLineLength + 2, fp) != NULL)
1690 lastTokenWasDash = NO;
1691 nNumbersOnThisLine = 0;
1693 if (GetToken (sumpToken, &tokenType, &sumpTokenP))
1695 /* printf ("%s (%d)\n", sumpToken, tokenType); */
1696 if (IsSame("[", sumpToken) == SAME)
1697 inSumpComment = YES;
1698 if (IsSame("]", sumpToken) == SAME)
1701 if (inSumpComment == NO)
1703 if (tokenType == NUMBER)
1705 sscanf (sumpToken, "%lf", &tempD);
1706 if (lastTokenWasDash == YES)
1708 nNumbersOnThisLine++;
1709 lastTokenWasDash = NO;
1711 else if (tokenType == DASH)
1713 lastTokenWasDash = YES;
1715 else if (tokenType != UNKNOWN_TOKEN_TYPE)
1718 lastTokenWasDash = NO;
1721 } while (*sumpToken);
1724 if (allDigitLine == NO)
1726 lastNonDigitLine = lineNum;
1731 if (nNumbersOnThisLine > 0)
1736 /* Now, check some aspects of the .p file. */
1737 if (inSumpComment == YES)
1739 MrBayesPrint ("%s Unterminated comment in file \"%s\"\n", spacer, fileName);
1742 if (numParamLines <= 0)
1744 MrBayesPrint ("%s No parameters were found in file or there characters not representing a number in last string of the file.\"%s\"\n", spacer, fileName);
1748 /* calculate burnin */
1749 if (chainParams.isSS == YES)
1755 if (chainParams.relativeBurnin == YES)
1756 burnin = (int) (chainParams.burninFraction * numParamLines);
1758 burnin = chainParams.chainBurnIn;
1761 /* check against burnin */
1762 if (burnin > numParamLines)
1764 MrBayesPrint ("%s No parameters can be sampled from file %s as the burnin (%d) exceeds the number of lines in last block (%d)\n",
1765 spacer, fileName, burnin, numParamLines);
1766 MrBayesPrint ("%s Try setting burnin to a number less than %d\n", spacer, numParamLines);
1770 /* Set some info in fileInfo */
1771 fileInfo->firstParamLine = lastNonDigitLine + burnin;
1772 fileInfo->headerLine = lastNonDigitLine;
1774 /* Calculate and check the number of columns and rows for the file; get header line at the same time */
1775 (void)fseek(fp, 0L, 0);
1776 for (lineNum=0; lineNum<lastNonDigitLine; lineNum++)
1777 if (fgets (s, fileInfo->longestLineLength + 2, fp)==NULL)
1779 strcpy(headerLine, s);
1780 for (; lineNum < lastNonDigitLine+burnin; lineNum++)
1781 if (fgets (s, fileInfo->longestLineLength + 2, fp)==NULL)
1787 while (fgets (s, fileInfo->longestLineLength + 2, fp) != NULL)
1791 lastTokenWasDash = NO;
1792 nNumbersOnThisLine = 0;
1794 if (GetToken (sumpToken, &tokenType, &sumpTokenP))
1796 if (IsSame("[", sumpToken) == SAME)
1797 inSumpComment = YES;
1798 if (IsSame("]", sumpToken) == SAME)
1800 if (inSumpComment == NO)
1802 if (tokenType == NUMBER)
1804 nNumbersOnThisLine++;
1805 lastTokenWasDash = NO;
1807 else if (tokenType == DASH)
1809 lastTokenWasDash = YES;
1811 else if (tokenType != UNKNOWN_TOKEN_TYPE)
1814 lastTokenWasDash = NO;
1817 } while (*sumpToken);
1819 if (allDigitLine == NO)
1821 MrBayesPrint ("%s Found a line with non-digit characters (line %d) in file %s\n", spacer, lineNum, fileName);
1826 if (nNumbersOnThisLine > 0)
1830 firstNumCols = nNumbersOnThisLine;
1833 if (nNumbersOnThisLine != firstNumCols)
1835 MrBayesPrint ("%s Number of columns is not even (%d in first line and %d in line %d of file %s)\n", spacer, firstNumCols, nNumbersOnThisLine, lineNum, fileName);
1842 fileInfo->numRows = nLines;
1843 fileInfo->numColumns = firstNumCols;
1845 /* set or check headers */
1846 if ((*headerNames) == NULL)
1848 GetHeaders (headerNames, headerLine, nHeaders);
1849 if (*nHeaders != fileInfo->numColumns)
1851 MrBayesPrint ("%s Expected %d headers but found %d headers\n", spacer, fileInfo->numColumns, *nHeaders);
1852 for (i=0; i<*nHeaders; i++)
1853 SafeFree ((void **) &((*headerNames)[i]));
1854 SafeFree ((void **) &(*headerNames));
1861 if (*nHeaders != fileInfo->numColumns)
1863 MrBayesPrint ("%s Expected %d columns but found %d columns\n", spacer, *nHeaders, fileInfo->numColumns);
1866 for (i=0, t=strtok(headerLine,"\t\n\r"); t!=NULL; t=strtok(NULL,"\t\n\r"), i++)
1870 MrBayesPrint ("%s Expected %d headers but found more headers.\n",
1871 spacer, fileInfo->numColumns);
1874 if (strcmp(t,(*headerNames)[i])!=0)
1876 MrBayesPrint ("%s Expected header '%s' for column %d but the header for this column was '%s' in file '%s'\n", spacer, (*headerNames)[i], i+1, t, fileName);
1877 MrBayesPrint ("%s It could be that some parameter values are not numbers and the whole string containing \n",spacer);
1878 MrBayesPrint ("%s this wrongly formated parameter is treated as a header.\n",spacer);
1884 MrBayesPrint ("%s Expected %d headers but found more headers.\n",spacer, fileInfo->numColumns);
1889 MrBayesPrint ("%s Expected header '%s' for column %d but the header for this column was '%s' in file '%s'\n",
1890 spacer, (*headerNames)[i], i+1, t, fileName);
1907 /***************************************************
1909 | FindHeader: Find token in list
1911 ----------------------------------------------------*/
1912 int FindHeader (char *token, char **headerNames, int nHeaders, int *index)
1914 int i, match=0, nMatches;
1918 for (i=0; i<nHeaders; i++)
1920 if (!strcmp(token,headerNames[i]))
1935 /* FreeParameterSamples: Free parameter samples space */
1936 void FreeParameterSamples (ParameterSample *parameterSamples)
1938 if (parameterSamples != NULL)
1940 free (parameterSamples[0].values[0]);
1941 free (parameterSamples[0].values);
1942 free (parameterSamples);
1947 /***************************************************
1949 | GetHeaders: Get headers from headerLine and put
1950 | them in list while updating nHeaders to reflect
1951 | the number of headers
1953 ----------------------------------------------------*/
1954 int GetHeaders (char ***headerNames, char *headerLine, int *nHeaders)
1959 for (s=strtok(headerLine," \t\n\r"); s!=NULL; s=strtok(NULL," \t\n\r"))
1961 if (AddString (headerNames, *nHeaders, s) == ERROR)
1963 MrBayesPrint ("%s Error adding header to list of headers \n", spacer, s);
1973 /* PrintMargLikes: Print marginal likelihoods to screen and to .lstat file */
1974 int PrintMargLikes (char *fileName, char **headerNames, int nHeaders, ParameterSample *parameterSamples, int nRuns, int nSamples)
1976 int i, j, len, longestHeader, *sampleCounts=NULL;
1981 /* calculate longest header */
1982 longestHeader = 9; /* length of 'parameter' */
1983 for (i=0; i<nHeaders; i++)
1985 strcpy (temp, headerNames[i]);
1986 len = (int) strlen(temp);
1987 for (j=0; modelIndicatorParams[j][0]!='\0'; j++)
1988 if (IsSame (temp,modelIndicatorParams[j]) != DIFFERENT)
1990 if (modelIndicatorParams[j][0]!='\0')
1992 if (!strcmp (temp, "Gen") || !strcmp (temp, "LnL") || !strcmp (temp, "LnPr"))
1994 if (len > longestHeader)
1995 longestHeader = len;
1998 /* open output file */
1999 strncpy (temp, fileName, 90);
2000 strcat (temp, ".pstat");
2001 fp = OpenNewMBPrintFile (temp);
2005 /* print unique identifier to the output file */
2006 if (strlen(stamp) > 1)
2007 MrBayesPrintf (fp, "[ID: %s]\n", stamp);
2009 /* allocate and set nSamples */
2010 sampleCounts = (int *) SafeCalloc (nRuns, sizeof(int));
2016 for (i=0; i<nRuns; i++)
2017 sampleCounts[i] = nSamples;
2019 /* print the header rows */
2020 MrBayesPrint ("\n");
2021 if (sumpParams.HPD == YES)
2022 MrBayesPrint ("%s %*c 95%% HPD Interval\n", spacer, longestHeader, ' ');
2024 MrBayesPrint ("%s %*c 95%% Cred. Interval\n", spacer, longestHeader, ' ');
2025 MrBayesPrint ("%s %*c --------------------\n", spacer, longestHeader, ' ');
2027 MrBayesPrint ("%s Parameter%*c Mean Variance Lower Upper Median", spacer, longestHeader-9, ' ');
2029 MrBayesPrint (" PSRF+ ");
2030 MrBayesPrint ("\n");
2032 MrBayesPrint ("%s ", spacer);
2033 for (j=0; j<longestHeader+1; j++)
2035 MrBayesPrint ("-----------------------------------------------------------");
2037 MrBayesPrint ("----------");
2038 MrBayesPrint ("\n");
2040 MrBayesPrintf (fp, "Parameter\tMean\tVariance\tLower\tUpper\tMedian\tPSRF\n");
2042 MrBayesPrintf (fp, "Parameter\tMean\tVariance\tLower\tUpper\tMedian\n");
2044 /* print table values */
2045 for (i=0; i<nHeaders; i++)
2047 strcpy (temp, headerNames[i]);
2048 len = (int) strlen(temp);
2049 for (j=0; modelIndicatorParams[j][0]!='\0'; j++)
2050 if (IsSame (temp,modelIndicatorParams[j]) != DIFFERENT)
2052 if (!strcmp (temp, "Gen") || !strcmp (temp, "LnL") || !strcmp (temp, "LnPr"))
2055 GetSummary (parameterSamples[i].values, nRuns, sampleCounts, &theStats, sumpParams.HPD);
2057 MrBayesPrint ("%s %-*s ", spacer, longestHeader, temp);
2058 MrBayesPrint ("%10.6lf %10.6lf %10.6lf %10.6lf %10.6lf", theStats.mean, theStats.var, theStats.lower, theStats.upper, theStats.median);
2059 MrBayesPrintf (fp, "%s", temp);
2060 MrBayesPrintf (fp, "\t%s", MbPrintNum(theStats.mean));
2061 MrBayesPrintf (fp, "\t%s", MbPrintNum(theStats.var));
2062 MrBayesPrintf (fp, "\t%s", MbPrintNum(theStats.lower));
2063 MrBayesPrintf (fp, "\t%s", MbPrintNum(theStats.upper));
2064 MrBayesPrintf (fp, "\t%s", MbPrintNum(theStats.median));
2067 if (theStats.PSRF < 0.0)
2069 MrBayesPrint (" NA ");
2070 MrBayesPrintf (fp, "NA");
2074 MrBayesPrint (" %7.3lf", theStats.PSRF);
2075 MrBayesPrintf (fp, "\t%s", MbPrintNum(theStats.PSRF));
2078 MrBayesPrint ("\n");
2079 MrBayesPrintf (fp, "\n");
2081 MrBayesPrint ("%s ", spacer);
2082 for (j=0; j<longestHeader+1; j++)
2084 MrBayesPrint ("-----------------------------------------------------------");
2086 MrBayesPrint ("----------");
2087 MrBayesPrint ("\n");
2090 MrBayesPrint ("%s + Convergence diagnostic (PSRF = Potential Scale Reduction Factor; Gelman\n", spacer);
2091 MrBayesPrint ("%s and Rubin, 1992) should approach 1.0 as runs converge.\n", spacer);
2095 free (sampleCounts);
2101 /* PrintModelStats: Print model stats to screen and to .mstat file */
2102 int PrintModelStats (char *fileName, char **headerNames, int nHeaders, ParameterSample *parameterSamples, int nRuns, int nSamples)
2104 int i, j, j1, j2, k, longestName, nElements, *modelCounts=NULL;
2105 MrBFlt f, *prob=NULL, *sum=NULL, *ssq=NULL, *min=NULL, *max=NULL, *stddev=NULL;
2108 ModelProb *elem = NULL;
2110 /* nHeaders - is a convenient synonym for number of column headers */
2112 /* check if we have any model indicator variables and also check for longest header */
2115 for (i=0; i<nHeaders; i++)
2117 for (j=0; strcmp(modelIndicatorParams[j],"")!=0; j++)
2119 if (IsSame (headerNames[i], modelIndicatorParams[j]) != DIFFERENT)
2122 for (j1=0; strcmp(modelElementNames[j][j1],"")!=0; j1++)
2124 j2 = (int)(strlen(headerNames[i]) + 2 + strlen(modelElementNames[j][j1]));
2125 if (j2 > longestName)
2133 /* return if nothing to do */
2137 /* open output file */
2138 MrBayesPrint ("%s Model probabilities above %1.3lf\n", spacer, sumpParams.minProb);
2139 MrBayesPrint ("%s Estimates saved to file \"%s.mstat\".\n", spacer, sumpParams.sumpOutfile);
2140 strncpy (temp,fileName,90);
2141 strcat (temp, ".mstat");
2142 fp = OpenNewMBPrintFile(temp);
2145 MrBayesPrint ("\n");
2147 /* print unique identifier to the output file */
2148 if (strlen(stamp) > 1)
2149 MrBayesPrintf (fp, "[ID: %s]\n", stamp);
2152 MrBayesPrintf (fp, "\n\n");
2155 MrBayesPrint ("%s %*c Posterior\n", spacer, longestName-5, ' ');
2156 MrBayesPrint ("%s Model%*c Probability\n", spacer, longestName-5, ' ');
2157 MrBayesPrint ("%s -----", spacer);
2158 for (i=0; i<longestName-5; i++)
2160 MrBayesPrint ("------------------\n");
2161 MrBayesPrintf (fp, "Model\tProbability\n");
2165 MrBayesPrint ("%s %*c Posterior Standard Min. Max. \n", spacer, longestName-5, ' ');
2166 MrBayesPrint ("%s Model%*c Probability Deviation Probability Probability\n", spacer, longestName-5, ' ');
2167 MrBayesPrint ("%s -----", spacer);
2168 for (i=0; i<longestName-5; i++)
2170 MrBayesPrint ("---------------------------------------------------------------\n");
2171 MrBayesPrintf (fp, "Model\tProbability\tStd_dev\tMin_prob\tMax_prob\n");
2174 /* calculate and print values */
2175 for (i=0; i<nHeaders; i++)
2177 for (j=0; modelIndicatorParams[j][0]!='\0'; j++)
2178 if (IsSame (headerNames[i], modelIndicatorParams[j]) != DIFFERENT)
2180 if (modelIndicatorParams[j][0] == '\0')
2183 for (nElements=0; modelElementNames[j][nElements][0]!='\0'; nElements++)
2186 modelCounts = (int *) SafeCalloc (nElements, sizeof(int));
2192 prob = (MrBFlt *) SafeCalloc (6*nElements, sizeof(MrBFlt));
2199 sum = prob + nElements;
2200 ssq = prob + 2*nElements;
2201 stddev = prob + 3*nElements;
2202 min = prob + 4*nElements;
2203 max = prob + 5*nElements;
2205 for (j1=0; j1<nElements; j1++)
2208 for (j1=0; j1<nRuns; j1++)
2210 for (j2=0; j2<nElements; j2++)
2211 modelCounts[j2] = 0;
2212 for (j2=0; j2<nSamples; j2++)
2213 modelCounts[(int)(parameterSamples[i].values[j1][j2] + 0.1)]++;
2214 for (j2=0; j2<nElements; j2++)
2216 f = (MrBFlt) modelCounts[j2] / (MrBFlt) nSamples;
2226 for (j1=0; j1<nElements; j1++)
2228 prob[j1] = sum[j1] / (MrBFlt) nRuns;
2229 f = ssq[j1] - (sum[j1] * sum[j1] / (MrBFlt) nRuns);
2234 stddev[j1] = sqrt (f);
2237 elem = (ModelProb *) SafeCalloc (nElements, sizeof(ModelProb));
2238 for (j1=0; j1<nElements; j1++)
2240 elem[j1].index = j1;
2241 elem[j1].prob = prob[j1];
2244 /* sort in terms of decreasing probabilities */
2245 qsort((void *)elem, (size_t)nElements, sizeof(ModelProb), CompareModelProbs);
2247 for (j1=0; j1<nElements; j1++)
2249 if (elem[j1].prob <= sumpParams.minProb)
2254 sprintf (temp, "%s[%s]", headerNames[i], modelElementNames[j][elem[j1].index]);
2255 MrBayesPrint ("%s %-*s %1.3lf\n", spacer, longestName, temp, prob[elem[j1].index]);
2256 MrBayesPrintf (fp, "%s\t%s\n", temp, MbPrintNum(prob[elem[j1].index]));
2258 else /* if (nRuns > 1) */
2260 sprintf (temp, "%s[%s]", headerNames[i], modelElementNames[j][elem[j1].index]);
2261 MrBayesPrint ("%s %-*s %1.3lf %1.3lf %1.3lf %1.3lf\n",
2262 spacer, longestName, temp, prob[elem[j1].index], stddev[elem[j1].index], min[elem[j1].index], max[elem[j1].index]);
2263 MrBayesPrintf (fp, "%s", temp);
2264 MrBayesPrintf (fp, "\t%s", MbPrintNum(prob[elem[j1].index]));
2265 MrBayesPrintf (fp, "\t%s", MbPrintNum(stddev[elem[j1].index]));
2266 MrBayesPrintf (fp, "\t%s", MbPrintNum(min[elem[j1].index]));
2267 MrBayesPrintf (fp, "\t%s", MbPrintNum(max[elem[j1].index]));
2268 MrBayesPrintf (fp, "\n");
2282 MrBayesPrint ("%s -----", spacer);
2283 for (i=0; i<longestName-5; i++)
2285 MrBayesPrint ("------------------\n\n");
2289 MrBayesPrint ("%s -----", spacer);
2290 for (i=0; i<longestName-5; i++)
2292 MrBayesPrint ("---------------------------------------------------------------\n\n");
2295 /* close output file */
2302 /* PrintOverlayPlot: Print overlay x-y plot of log likelihood vs. generation for several runs */
2303 int PrintOverlayPlot (MrBFlt **xVals, MrBFlt **yVals, int nRuns, int startingFrom, int nSamples)
2305 int i, j, k, k2, n, screenHeight, screenWidth, numY[60], width;
2306 char plotSymbol[15][60];
2307 MrBFlt x, y, minX, maxX, minY, maxY, meanY[60];
2310 MrBayesPrint ("\n%s Overlay plot for both runs:\n", spacer);
2312 MrBayesPrint ("\n%s Overlay plot for all %d runs:\n", spacer, sumpParams.numRuns);
2314 MrBayesPrint ("%s (1 = Run number 1; 2 = Run number 2 etc.; x = Run number 10 or above; * = Several runs)\n", spacer);
2316 MrBayesPrint ("%s (1 = Run number 1; 2 = Run number 2 etc.; * = Several runs)\n", spacer);
2318 MrBayesPrint ("%s (1 = Run number 1; 2 = Run number 2; * = Both runs)\n", spacer);
2320 /* print overlay x-y plot of log likelihood vs. generation for all runs */
2321 screenWidth = 60; /* don't change this without changing numY, meanY, and plotSymbol declared above */
2324 /* find minX, minY, maxX, and maxY over all runs */
2325 minX = minY = 1000000000.0;
2326 maxX = maxY = -1000000000.0;
2327 for (n=0; n<nRuns; n++)
2329 for (i=startingFrom; i<startingFrom+nSamples; i++)
2338 for (n=0; n<nRuns; n++)
2343 for (i=startingFrom; i<startingFrom+nSamples; i++)
2346 k = (int) (((x - minX) / (maxX - minX)) * screenWidth);
2349 if (k >= screenWidth)
2350 k = screenWidth - 1;
2378 /* initialize the plot symbols */
2379 for (i=0; i<screenHeight; i++)
2380 for (j=0; j<screenWidth; j++)
2381 plotSymbol[i][j] = ' ';
2383 /* assemble the plot symbols */
2384 for (n=0; n<nRuns; n++)
2386 /* find the plot points for this run */
2387 for (i=0; i<screenWidth; i++)
2392 for (i=startingFrom; i<startingFrom+nSamples; i++)
2396 k = (int)(((x - minX) / (maxX - minX)) * screenWidth);
2397 if (k >= screenWidth)
2398 k = screenWidth - 1;
2405 /* transfer these points to the overlay */
2406 for (i=0; i<screenWidth; i++)
2410 k = (int) ((((meanY[i] / numY[i]) - minY)/ (maxY - minY)) * screenHeight);
2413 else if (k >= screenHeight)
2414 k = screenHeight - 1;
2415 if (plotSymbol[k][i] == ' ')
2418 plotSymbol[k][i] = '1' + n;
2420 plotSymbol[k][i] = 'x';
2423 plotSymbol[k][i] = '*';
2428 /* now print the overlay plot */
2429 MrBayesPrint ("\n%s +", spacer);
2430 for (i=0; i<screenWidth; i++)
2432 MrBayesPrint ("+ %1.2lf\n", maxY);
2433 for (i=screenHeight-1; i>=0; i--)
2435 MrBayesPrint ("%s |", spacer);
2436 for (j=0; j<screenWidth; j++)
2438 MrBayesPrint ("%c", plotSymbol[i][j]);
2440 MrBayesPrint ("|\n");
2442 MrBayesPrint ("%s +", spacer);
2443 for (i=0; i<screenWidth; i++)
2445 if (i % (screenWidth/10) == 0 && i != 0)
2450 MrBayesPrint ("+ %1.2lf\n", minY);
2451 MrBayesPrint ("%s ^", spacer);
2452 for (i=0; i<screenWidth; i++)
2454 MrBayesPrint ("^\n");
2455 MrBayesPrint ("%s %1.0lf", spacer, minX);
2457 width=(int)(log10(minX));
2458 else if ((int)minX==0)
2461 width=(int)(log10(-minX))+1;
2462 for (i=0; i<screenWidth-width; i++)
2464 MrBayesPrint ("%1.0lf\n\n", maxX);
2470 /* PrintParamStats: Print parameter table (not model indicator params) to screen and .pstat file */
2471 int PrintParamStats (char *fileName, char **headerNames, int nHeaders, ParameterSample *parameterSamples, int nRuns, int nSamples)
2473 int i, j, len, longestHeader, *sampleCounts=NULL;
2474 static char *temp=NULL;
2479 /* calculate longest header */
2480 longestHeader = 9; /* length of 'parameter' */
2481 for (i=0; i<nHeaders; i++)
2483 SafeStrcpy (&temp, headerNames[i]);
2484 len = (int) strlen(temp);
2485 for (j=0; modelIndicatorParams[j][0]!='\0'; j++)
2486 if (IsSame (temp,modelIndicatorParams[j]) != DIFFERENT)
2488 if (modelIndicatorParams[j][0]!='\0')
2490 if (!strcmp (temp, "Gen") || !strcmp (temp, "LnL") || !strcmp (temp, "LnPr"))
2492 if (len > longestHeader)
2493 longestHeader = len;
2496 /* open output file */
2497 strncpy (tempf, fileName, 90);
2498 strcat (tempf, ".pstat");
2499 fp = OpenNewMBPrintFile (tempf);
2503 /* print unique identifier to the output file */
2504 if (strlen(stamp) > 1)
2505 MrBayesPrintf (fp, "[ID: %s]\n", stamp);
2507 /* allocate and set nSamples */
2508 sampleCounts = (int *) SafeCalloc (nRuns, sizeof(int));
2514 for (i=0; i<nRuns; i++)
2515 sampleCounts[i] = nSamples;
2517 /* print the header rows */
2518 MrBayesPrint ("\n");
2519 if (sumpParams.HPD == YES)
2520 MrBayesPrint ("%s %*c 95%% HPD Interval\n", spacer, longestHeader, ' ');
2522 MrBayesPrint ("%s %*c 95%% Cred. Interval\n", spacer, longestHeader, ' ');
2523 MrBayesPrint ("%s %*c --------------------\n", spacer, longestHeader, ' ');
2526 MrBayesPrint ("%s Parameter%*c Mean Variance Lower Upper Median min ESS* avg ESS PSRF+ ", spacer, longestHeader-9, ' ');
2528 MrBayesPrint ("%s Parameter%*c Mean Variance Lower Upper Median ESS*", spacer, longestHeader-9, ' ');
2529 MrBayesPrint ("\n");
2531 MrBayesPrint ("%s ", spacer);
2532 for (j=0; j<longestHeader+1; j++)
2534 MrBayesPrint ("---------------------------------------------------------------------");
2536 MrBayesPrint ("-------------------");
2537 MrBayesPrint ("\n");
2539 MrBayesPrintf (fp, "Parameter\tMean\tVariance\tLower\tUpper\tMedian\tminESS\tavgESS\tPSRF\n");
2541 MrBayesPrintf (fp, "Parameter\tMean\tVariance\tLower\tUpper\tMedian\tESS\n");
2543 /* print table values */
2544 for (i=0; i<nHeaders; i++)
2546 SafeStrcpy(&temp, headerNames[i]);
2547 len = (int) strlen(temp);
2548 for (j=0; modelIndicatorParams[j][0]!='\0'; j++)
2549 if (IsSame (temp,modelIndicatorParams[j]) != DIFFERENT)
2551 if (modelIndicatorParams[j][0]!='\0')
2553 if (!strcmp (temp, "Gen") || !strcmp (temp, "LnL") || !strcmp (temp, "LnPr"))
2556 GetSummary (parameterSamples[i].values, nRuns, sampleCounts, &theStats, sumpParams.HPD);
2558 MrBayesPrint ("%s %-*s ", spacer, longestHeader, temp);
2559 MrBayesPrint ("%10.6lf %10.6lf %10.6lf %10.6lf %10.6lf", theStats.mean, theStats.var, theStats.lower, theStats.upper, theStats.median);
2560 MrBayesPrintf (fp, "%s", temp);
2561 MrBayesPrintf (fp, "\t%s", MbPrintNum(theStats.mean));
2562 MrBayesPrintf (fp, "\t%s", MbPrintNum(theStats.var));
2563 MrBayesPrintf (fp, "\t%s", MbPrintNum(theStats.lower));
2564 MrBayesPrintf (fp, "\t%s", MbPrintNum(theStats.upper));
2565 MrBayesPrintf (fp, "\t%s", MbPrintNum(theStats.median));
2567 if (theStats.minESS == theStats.minESS)
2569 MrBayesPrintf (fp, "\t%s", MbPrintNum(theStats.minESS));
2570 MrBayesPrint (" %8.2lf", theStats.minESS);
2574 MrBayesPrint (" NA ");
2575 MrBayesPrintf (fp, "NA");
2579 if (theStats.minESS == theStats.minESS)
2581 MrBayesPrint (" %8.2lf", theStats.avrESS);
2582 MrBayesPrintf (fp, "\t%s", MbPrintNum(theStats.avrESS));
2586 MrBayesPrint (" NA ");
2587 MrBayesPrintf (fp, "NA");
2589 if (theStats.PSRF < 0.0)
2591 MrBayesPrint (" NA ");
2592 MrBayesPrintf (fp, "NA");
2596 MrBayesPrint (" %7.3lf", theStats.PSRF);
2597 MrBayesPrintf (fp, "\t%s", MbPrintNum(theStats.PSRF));
2600 MrBayesPrint ("\n");
2601 MrBayesPrintf (fp, "\n");
2603 MrBayesPrint ("%s ", spacer);
2604 for (j=0; j<longestHeader+1; j++)
2606 MrBayesPrint ("---------------------------------------------------------------------");
2608 MrBayesPrint ("-------------------");
2609 MrBayesPrint ("\n");
2612 MrBayesPrint ("%s * Convergence diagnostic (ESS = Estimated Sample Size); min and avg values\n", spacer);
2613 MrBayesPrint ("%s correspond to minimal and average ESS among runs. \n", spacer);
2614 MrBayesPrint ("%s ESS value below 100 may indicate that the parameter is undersampled. \n", spacer);
2615 MrBayesPrint ("%s + Convergence diagnostic (PSRF = Potential Scale Reduction Factor; Gelman\n", spacer);
2616 MrBayesPrint ("%s and Rubin, 1992) should approach 1.0 as runs converge.\n", spacer);
2620 MrBayesPrint ("%s * Convergence diagnostic (ESS = Estimated Sample Size); ESS value \n", spacer);
2621 MrBayesPrint ("%s below 100 may indicate that the parameter is undersampled. \n", spacer);
2623 MrBayesPrint ("\n\n");
2626 free (sampleCounts);
2627 SafeFree ((void **)&temp);
2633 /* PrintPlot: Print x-y plot of log likelihood vs. generation */
2634 int PrintPlot (MrBFlt *xVals, MrBFlt *yVals, int numVals)
2636 int i, j, k, numY[60], screenWidth, screenHeight;
2637 MrBFlt x, y, minX, maxX, minY, maxY, meanY[60], diff;
2639 /* print x-y plot of log likelihood vs. generation */
2640 screenWidth = 60; /* don't change this without changing numY and meanY, declared above */
2643 /* find minX and maxX */
2646 for (i=0; i<numVals; i++)
2655 /* collect Y data */
2656 for (i=0; i<screenWidth; i++)
2661 for (i=0; i<numVals; i++)
2665 k = (int)(((x - minX) / (maxX - minX)) * screenWidth);
2666 if (k >= screenWidth)
2667 k = screenWidth - 1;
2674 /* find minY and maxY */
2675 minY = maxY = meanY[0] / numY[0];
2676 for (i=0; i<screenWidth; i++)
2678 if (meanY[i] == 0) /* with some compilers if (NaN < 1) is equal true !!! so we realy need this check*/
2680 meanY[i] /= numY[i];
2681 if (meanY[i] < minY)
2683 if (meanY[i] > maxY)
2687 /* find difference */
2691 MrBayesPrint ("\n +");
2692 for (i=0; i<screenWidth; i++)
2694 MrBayesPrint ("+ %1.3lf\n", maxY);
2695 for (j=screenHeight-1; j>=0; j--)
2697 MrBayesPrint (" |");
2698 for (i=0; i<screenWidth; i++)
2702 if ((meanY[i] > ((diff/screenHeight)*j)+minY && meanY[i] <= ((diff/screenHeight)*(j+1))+minY) ||
2703 (j == 0 && meanY[i] <= minY))
2713 MrBayesPrint ("|\n");
2715 MrBayesPrint (" +");
2716 for (i=0; i<screenWidth; i++)
2718 if (i % (screenWidth/10) == 0 && i != 0)
2723 MrBayesPrint ("+ %1.3lf\n", minY);
2724 MrBayesPrint (" ^");
2725 for (i=0; i<screenWidth; i++)
2727 MrBayesPrint ("^\n");
2728 MrBayesPrint (" %1.0lf", minX);
2732 j = (int)(log10(minX)) + 1;
2733 for (i=0; i<screenWidth-j; i++)
2735 MrBayesPrint ("%1.0lf\n\n", maxX);
2741 void PrintPlotHeader (void)
2743 MrBayesPrint ("\n");
2744 if (sumpParams.numRuns > 1)
2746 MrBayesPrint ("%s Below are rough plots of the generation (x-axis) versus the log \n", spacer);
2747 MrBayesPrint ("%s probability of observing the data (y-axis). You can use these \n", spacer);
2748 MrBayesPrint ("%s graphs to determine what the burn in for your analysis should be. \n", spacer);
2749 MrBayesPrint ("%s When the log probability starts to plateau you may be at station- \n", spacer);
2750 MrBayesPrint ("%s arity. Sample trees and parameters after the log probability \n", spacer);
2751 MrBayesPrint ("%s plateaus. Of course, this is not a guarantee that you are at sta- \n", spacer);
2752 MrBayesPrint ("%s tionarity. Also examine the convergence diagnostics provided by \n", spacer);
2753 MrBayesPrint ("%s the 'sump' and 'sumt' commands for all the parameters in your \n", spacer);
2754 MrBayesPrint ("%s model. Remember that the burn in is the number of samples to dis- \n", spacer);
2755 MrBayesPrint ("%s card. There are a total of ngen / samplefreq samples taken during \n", spacer);
2756 MrBayesPrint ("%s a MCMC analysis. \n", spacer);
2760 MrBayesPrint ("%s Below is a rough plot of the generation (x-axis) versus the log \n", spacer);
2761 MrBayesPrint ("%s probability of observing the data (y-axis). You can use this \n", spacer);
2762 MrBayesPrint ("%s graph to determine what the burn in for your analysis should be. \n", spacer);
2763 MrBayesPrint ("%s When the log probability starts to plateau you may be at station- \n", spacer);
2764 MrBayesPrint ("%s arity. Sample trees and parameters after the log probability \n", spacer);
2765 MrBayesPrint ("%s plateaus. Of course, this is not a guarantee that you are at sta- \n", spacer);
2766 MrBayesPrint ("%s analysis should be. When the log probability starts to plateau \n", spacer);
2767 MrBayesPrint ("%s tionarity. When possible, run multiple analyses starting from dif-\n", spacer);
2768 MrBayesPrint ("%s ferent random trees; if the inferences you make for independent \n", spacer);
2769 MrBayesPrint ("%s analyses are the same, this is reasonable evidence that the chains\n", spacer);
2770 MrBayesPrint ("%s have converged. You can use MrBayes to run several independent \n", spacer);
2771 MrBayesPrint ("%s analyses simultaneously. During such a run, MrBayes will monitor \n", spacer);
2772 MrBayesPrint ("%s the convergence of topologies. After the run has been completed, \n", spacer);
2773 MrBayesPrint ("%s the 'sumt' and 'sump' functions will provide additional conver- \n", spacer);
2774 MrBayesPrint ("%s gence diagnostics for all the parameters in your model. Remember \n", spacer);
2775 MrBayesPrint ("%s that the burn in is the number of samples to discard. There are \n", spacer);
2776 MrBayesPrint ("%s a total of ngen / samplefreq samples taken during a MCMC analysis.\n", spacer);
2781 /* ReadParamSamples: Read parameter samples from .p file */
2782 int ReadParamSamples (char *fileName, SumpFileInfo *fileInfo, ParameterSample *parameterSamples, int runNo)
2784 char sumpToken[CMD_STRING_LENGTH], *s=NULL, *p;
2785 int inSumpComment, lineNum, numLinesRead, numLinesToRead, column, lastTokenWasDash,
2791 if ((fp = OpenTextFileR (fileName)) == NULL)
2794 /* allocate space for reading lines */
2795 s = (char *) SafeCalloc (fileInfo->longestLineLength + 10, sizeof(char));
2797 /* fast forward to beginning of last unburned parameter line. */
2798 for (lineNum=0; lineNum<fileInfo->firstParamLine; lineNum++)
2799 if (fgets (s, fileInfo->longestLineLength + 5, fp) == 0)
2802 /* parse file, line-by-line. We are only parsing lines that have digits that should be read. */
2804 numLinesToRead = fileInfo->numRows;
2806 while (fgets (s, fileInfo->longestLineLength + 1, fp) != NULL)
2808 lastTokenWasDash = NO;
2812 if (GetToken (sumpToken, &tokenType, &p))
2814 if (IsSame("[", sumpToken) == SAME)
2815 inSumpComment = YES;
2816 if (IsSame("]", sumpToken) == SAME)
2818 if (inSumpComment == NO)
2820 if (tokenType == NUMBER)
2822 /* read the number */
2823 if (column >= fileInfo->numColumns)
2825 MrBayesPrint ("%s Too many values read on line %d of file %s\n", spacer, lineNum, fileName);
2828 sscanf (sumpToken, "%lf", &tempD);
2829 if (lastTokenWasDash == YES)
2831 parameterSamples[column].values[runNo][numLinesRead] = tempD;
2833 lastTokenWasDash = NO;
2835 else if (tokenType == DASH)
2837 lastTokenWasDash = YES;
2839 else if (tokenType != UNKNOWN_TOKEN_TYPE)
2841 /* we have a problem */
2842 MrBayesPrint ("%s Line %d of file %s has non-digit characters\n", spacer, lineNum, fileName);
2846 } while (*sumpToken);
2849 if (column == fileInfo->numColumns)
2851 else if (column != 0)
2853 MrBayesPrint ("%s Too few values on line %d of file %s\n", spacer, lineNum, fileName);
2859 /* Check how many parameter line was read in. */
2860 if (numLinesRead != numLinesToRead)
2862 MrBayesPrint ("%s Unable to read all lines that should contain parameter samples\n", spacer);
2880 /* the following are moved from sumt.c to combine with sump.c */
2881 PartCtr *AddSumtPartition (PartCtr *r, PolyTree *t, PolyNode *p, int runId)
2883 int i, n, comp, nLongsNeeded = sumtParams.BitsLongsNeeded;
2888 /* create a new node */
2889 r = AllocPartCtr ();
2892 numUniqueSplitsFound++;
2893 for (i=0; i<nLongsNeeded; i++)
2894 r->partition[i] = p->partition[i];
2895 for (i=0; i<sumtParams.numRuns; i++)
2897 r->left = r->right = NULL;
2899 if (sumtParams.brlensDef == YES)
2900 r->length[runId][0]= p->length;
2901 if (sumtParams.isClock == YES)
2902 r->height[runId][0]= p->depth;
2903 if (sumtParams.isCalibrated == YES)
2904 r->age[runId][0]= p->age;
2905 for (i=0; i<sumtParams.nESets; i++)
2906 r->nEvents[i][runId][0] = t->nEvents[i][p->index];
2907 for (i=0; i<sumtParams.nBSets; i++)
2909 r->bLen [i][runId][0] = t->effectiveBrLen[i][p->index];
2910 r->bRate[i][runId][0] = t->effectiveBrLen[i][p->index] / p->length;
2912 if (t->popSizeSet == YES)
2913 r->popSize[runId][0] = t->popSize[p->index];
2919 for (i=0; i<nLongsNeeded; i++)
2921 if (r->partition[i] != p->partition[i])
2925 if (i == nLongsNeeded)
2927 else if (r->partition[i] < p->partition[i])
2932 if (comp == 0) /* repeated partition */
2934 n = r->count[runId];
2935 /* check if we need to allocate more space */
2936 if (n % ALLOC_LEN == 0)
2938 /* allocate more space */
2939 if (sumtParams.brlensDef == YES)
2940 r->length[runId] = (MrBFlt *) SafeRealloc ((void *)r->length[runId], ((size_t)n+ALLOC_LEN)*sizeof(MrBFlt));
2941 if (sumtParams.isClock == YES)
2942 r->height[runId] = (MrBFlt *) SafeRealloc ((void *)r->height[runId], ((size_t)n+ALLOC_LEN)*sizeof(MrBFlt));
2943 if (sumtParams.isCalibrated == YES)
2944 r->age[runId] = (MrBFlt *) SafeRealloc ((void *)r->age[runId], ((size_t)n+ALLOC_LEN)*sizeof(MrBFlt));
2945 if (sumtParams.nESets > 0)
2947 for (i=0; i<sumtParams.nESets; i++)
2948 r->nEvents[i][runId] = (int *) SafeRealloc ((void *)r->nEvents[i][runId], ((size_t)n+ALLOC_LEN)*sizeof(int));
2950 if (sumtParams.nBSets > 0)
2952 for (i=0; i<sumtParams.nBSets; i++)
2954 r->bRate[i][runId] = (MrBFlt *) SafeRealloc ((void *)r->bRate[i][runId], ((size_t)n+ALLOC_LEN)*sizeof(MrBFlt));
2955 r->bLen [i][runId] = (MrBFlt *) SafeRealloc ((void *)r->bLen [i][runId], ((size_t)n+ALLOC_LEN)*sizeof(MrBFlt));
2958 if (sumtParams.popSizeSet == YES)
2959 r->popSize[runId] = (MrBFlt *) SafeRealloc ((void *)r->popSize[runId], ((size_t)n+ALLOC_LEN)*sizeof(MrBFlt));
2964 if (sumtParams.brlensDef == YES)
2965 r->length[runId][n]= p->length;
2966 if (sumtParams.isClock == YES)
2967 r->height[runId][n]= p->depth;
2968 if (sumtParams.isCalibrated == YES)
2969 r->age[runId][n]= p->age;
2970 if (sumtParams.nESets > 0)
2972 for (i=0; i<sumtParams.nESets; i++)
2973 r->nEvents[i][runId][n] = t->nEvents[i][p->index];
2975 if (sumtParams.nBSets > 0)
2977 for (i=0; i<sumtParams.nBSets; i++)
2979 r->bLen [i][runId][n] = t->effectiveBrLen[i][p->index];
2980 r->bRate[i][runId][n] = t->effectiveBrLen[i][p->index] / p->length;
2983 if (sumtParams.popSizeSet == YES)
2984 r->popSize[runId][n] = t->popSize[p->index];
2986 else if (comp < 0) /* greater than -> into left subtree */
2988 if ((r->left = AddSumtPartition (r->left, t, p, runId)) == NULL)
2996 /* smaller than -> into right subtree */
2997 if ((r->right = AddSumtPartition (r->right, t, p, runId)) == NULL)
3009 TreeCtr *AddSumtTree (TreeCtr *r, int *order)
3016 /* create a new node */
3020 numUniqueTreesFound++;
3021 for (i=0; i<sumtParams.orderLen; i++)
3022 r->order[i] = order[i];
3027 for (i=0; i<sumtParams.orderLen; i++)
3028 if (r->order[i] != order[i])
3031 if (i==sumtParams.orderLen)
3033 else if (order[i] < r->order[i])
3038 if (comp == 0) /* repeated partition */
3040 else if (comp < 0) /* greater than -> into left subtree */
3042 if ((r->left = AddSumtTree (r->left, order)) == NULL)
3050 /* smaller than -> into right subtree */
3051 if ((r->right = AddSumtTree (r->right, order)) == NULL)
3063 /* AllocPartCtr: Allocate space for one partition counter node using info in sumtParams */
3064 PartCtr *AllocPartCtr ()
3069 /* allocate basic stuff */
3070 r = (PartCtr *) SafeCalloc (1, sizeof(PartCtr));
3071 r->left = r->right = NULL;
3072 r->partition = (BitsLong *) SafeCalloc ((size_t)(sumtParams.BitsLongsNeeded), sizeof(BitsLong));
3073 r->count = (int *) SafeCalloc ((size_t)(sumtParams.numRuns), sizeof (int));
3074 if (sumtParams.brlensDef)
3076 r->length = (MrBFlt **) SafeCalloc ((size_t)(sumtParams.numRuns), sizeof (MrBFlt *));
3077 for (i=0; i<sumtParams.numRuns; i++)
3078 r->length[i] = (MrBFlt *) SafeCalloc ((size_t)ALLOC_LEN, sizeof(MrBFlt));
3080 if (sumtParams.isClock)
3082 r->height = (MrBFlt **) SafeCalloc ((size_t)(sumtParams.numRuns), sizeof (MrBFlt *));
3083 for (i=0; i<sumtParams.numRuns; i++)
3084 r->height[i] = (MrBFlt *) SafeCalloc ((size_t)ALLOC_LEN, sizeof(MrBFlt));
3085 r->age = (MrBFlt **) SafeCalloc ((size_t)(sumtParams.numRuns), sizeof (MrBFlt *));
3086 for (i=0; i<sumtParams.numRuns; i++)
3087 r->age[i] = (MrBFlt *) SafeCalloc ((size_t)ALLOC_LEN, sizeof(MrBFlt));
3090 /* allocate relaxed clock parameters: eRate, nEvents, bRate */
3091 if (sumtParams.nESets > 0)
3092 r->nEvents = (int ***) SafeCalloc ((size_t)(sumtParams.nESets), sizeof(int **));
3093 for (i=0; i<sumtParams.nESets; i++)
3095 r->nEvents[i] = (int **) SafeCalloc ((size_t)(sumtParams.numRuns), sizeof(int *));
3096 for (j=0; j<sumtParams.numRuns; j++)
3097 r->nEvents[i][j] = (int *) SafeCalloc ((size_t) ALLOC_LEN, sizeof(int));
3099 if (sumtParams.nBSets > 0)
3101 r->bLen = (MrBFlt ***) SafeCalloc ((size_t)(sumtParams.nBSets), sizeof(MrBFlt **));
3102 r->bRate = (MrBFlt ***) SafeCalloc ((size_t)(sumtParams.nBSets), sizeof(MrBFlt **));
3104 for (i=0; i<sumtParams.nBSets; i++)
3106 r->bLen[i] = (MrBFlt **) SafeCalloc ((size_t)(sumtParams.numRuns), sizeof(MrBFlt *));
3107 r->bRate[i] = (MrBFlt **) SafeCalloc ((size_t)(sumtParams.numRuns), sizeof(MrBFlt *));
3108 for (j=0; j<sumtParams.numRuns; j++)
3110 r->bLen[i][j] = (MrBFlt *) SafeCalloc ((size_t)ALLOC_LEN, sizeof(MrBFlt));
3111 r->bRate[i][j] = (MrBFlt *) SafeCalloc ((size_t)ALLOC_LEN, sizeof(MrBFlt));
3114 if (sumtParams.popSizeSet == YES)
3116 r->popSize = (MrBFlt **) SafeCalloc ((size_t)(sumtParams.numRuns), sizeof (MrBFlt *));
3117 for (i=0; i<sumtParams.numRuns; i++)
3118 r->popSize[i] = (MrBFlt *) SafeCalloc ((size_t)ALLOC_LEN, sizeof(MrBFlt));
3125 /* AllocTreeCtr: Allocate space for a tree counter node using info in sumtParams struct*/
3126 TreeCtr *AllocTreeCtr ()
3130 r = (TreeCtr *) SafeCalloc (1, sizeof(TreeCtr));
3132 r->left = r->right = NULL;
3134 r->order = (int *) SafeCalloc ((size_t)(sumtParams.orderLen), sizeof(int));
3140 void CalculateTreeToTreeDistance (Tree *tree1, Tree *tree2, MrBFlt *d1, MrBFlt *d2, MrBFlt *d3)
3143 MrBFlt treeLen1=0.0, treeLen2=0.0;
3144 TreeNode *p, *q=NULL;
3146 (*d1) = (*d2) = (*d3) = 0.0;
3148 /* set distance-based measures to max value */
3149 if (sumtParams.brlensDef == YES)
3151 treeLen1 = TreeLen(tree1);
3152 treeLen2 = TreeLen(tree2);
3153 (*d2) = treeLen1 + treeLen2;
3157 /* now we can get distances in a single pass */
3158 for (i=0; i<tree1->nNodes; i++)
3160 p = tree1->allDownPass[i];
3161 for (j=0; j<tree2->nNodes; j++)
3163 q = tree2->allDownPass[j];
3164 for (k=0; k<sumtParams.BitsLongsNeeded; k++)
3165 if (p->partition[k] != q->partition[k])
3167 if (k == sumtParams.BitsLongsNeeded)
3170 if (j < tree2->nNodes)
3173 if (sumtParams.brlensDef == YES)
3175 (*d2) -= (p->length + q->length - fabs(p->length - q->length));
3176 (*d3) -= (p->length/treeLen1 + q->length/treeLen2 - fabs(p->length/treeLen1 - q->length/treeLen2));
3179 else /* if (k < sumtParams.BitsLongsNeeded) */
3187 printf ("DISTANCES: %lf %lf %lf (%lf %lf)\n", *d1, *d2, *d3, tl1, tl2);
3188 for (i=0; i<nnds; i++)
3190 printf ("%4d -- %4d (%lf) %4d (%lf)\n", i, list1[i], lengths1[i], list2[i], lengths2[i]);
3196 /* ConTree: Construct consensus tree */
3197 int ConTree (PartCtr **treeParts, int numTreeParts)
3199 int i, j, targetNode, nBits, isCompat, numTerminalsEncountered;
3200 BitsLong x, *partition = NULL, bitsLongOne;
3201 MrBFlt freq, freqInterapted=0;
3202 PolyTree *t, *t2=NULL;
3203 PolyNode *p, *q, *r, *ql, *pl;
3206 int isFirstLoop=1, isInterapted=0;
3208 /* we use a BitsLong variable set to one to escape dependence on interpretation of constant (signed?) int/long '1L' */
3211 /* check that we have at least three species */
3212 if (sumtParams.numTaxa < 3)
3214 MrBayesPrint ("%s Too few taxa included to show consensus trees\n", spacer);
3219 /* now, make a consensus tree */
3220 /* first allocate and initialize consensus tree */
3221 t = AllocatePolyTree(sumtParams.numTaxa);
3224 MrBayesPrint ("%s Could not allocate consensus tree\n", spacer);
3227 t->isRooted = sumtParams.isRooted;
3228 t->isClock = sumtParams.isClock;
3229 t->isRelaxed = sumtParams.isRelaxed;
3231 /* initialize consensus tree nodes */
3232 for (i=0; i<sumtParams.numTaxa; i++)
3234 t->nodes[i].left = NULL;
3235 t->nodes[i].sib = NULL;
3236 t->nodes[i].index = i;
3237 t->nodes[i].partitionIndex = -1; /* partition ID */
3238 t->nodes[i].age = 0.0; /* temporally set to minimum value to allow any insertion in front of the terminal before actual
3239 values of age and depth are available */
3240 strcpy(t->nodes[i].label, sumtParams.taxaNames[i]);
3241 t->nodes[i].depth = 0.0;
3243 for (; i<t->memNodes; i++)
3245 t->nodes[i].left = NULL;
3246 t->nodes[i].sib = NULL;
3247 t->nodes[i].index = i;
3248 t->nodes[i].partitionIndex = -1; /* partition ID */
3249 strcpy (t->nodes[i].label, "");
3253 ->x counts number of subtended terminals
3254 make sure t->root->left is in outgroup */
3255 p = t->root = &t->nodes[sumtParams.numTaxa];
3256 p->anc = p->sib = NULL;
3257 p->x = sumtParams.numTaxa;
3258 p->age = MRBFLT_MAX; /* temporarily set to maximum value to allow any insertion in front of the root before actual values of age and depth are available */
3259 p->depth = MRBFLT_MAX;
3265 for (i=0; i<sumtParams.numTaxa; i++)
3269 q->sib = &t->nodes[i];
3277 /* Resolve bush according to partitions.
3278 Partitions may include incompatible ones.
3279 Partitions must be sorted from most frequent to least frequent
3280 for quit test to work when a 50% majority rule tree is requested
3281 and in general for consensus tree to be correct. */
3282 t->nNodes = sumtParams.numTaxa + 1;
3284 if (sumtParams.isRooted == YES)
3285 targetNode = 2 * sumtParams.numTaxa - 2;
3287 targetNode = 2 * sumtParams.numTaxa - 3;
3289 numTerminalsEncountered = 0;
3290 for (i=0; i<numUniqueSplitsFound; i++)
3293 part = treeParts[i];
3295 /* calculate frequency and test if time to quit */
3296 if (t->nNodes > targetNode && numTerminalsEncountered == sumtParams.numTaxa)
3298 freq = (MrBFlt)(part->totCount) / (MrBFlt)(sumtParams.numTreesSampled);
3299 if (freq < 0.50 && !strcmp(sumtParams.sumtConType, "Halfcompat"))
3303 partition = part->partition;
3305 /* count bits in this partition */
3306 for (j=nBits=0; j<sumtParams.BitsLongsNeeded; j++)
3309 for (x = partition[j]; x != 0; x &= (x - 1))
3313 /* find out if this is an informative partition */
3314 if (nBits == sumtParams.numTaxa || nBits == 0)
3316 /* this is the root (for setting age of root node when tree is dated) */
3318 q->partitionIndex = i;
3319 if (sumtParams.isClock == YES)
3321 GetSummary(part->height, sumtParams.numRuns, part->count, &theStats, sumtParams.HPD);
3322 q->depth = theStats.median;
3323 for (p = q->left; p!=NULL; p = p->sib)
3325 if (q->depth <= p->depth)
3328 assert (p==NULL); /* Root always has 100% freq and it should be older than any other node that has 100% freq. */
3330 if (sumtParams.isCalibrated == YES)
3332 GetSummary(part->age, sumtParams.numRuns, part->count, &theStats, sumtParams.HPD);
3333 q->age = theStats.median;
3334 for (p = q->left; p!=NULL; p = p->sib)
3336 if (q->age <= p->age)
3339 assert (p==NULL); /* Root always has 100% freq and it should be older than any other node that has 100% freq. */
3342 else if (nBits > 1 && !(nBits == sumtParams.numTaxa - 1 && sumtParams.isRooted == NO))
3344 /* this is an informative partition */
3345 /* find anc of partition */
3346 j = FirstTaxonInPartition (partition, sumtParams.BitsLongsNeeded);
3347 for (p = &t->nodes[j]; p!=NULL; p = p->anc)
3351 /* do not include if incompatible with ancestor or any of descendants
3352 do not check terminals or root because it is
3353 redundant and partitions have not necessarily been set for those */
3355 if (p->anc != NULL && IsPartNested(partition, p->partition, sumtParams.BitsLongsNeeded)==NO)
3359 for (q=p->left; q!=NULL; q=q->sib)
3361 if (q->x > 1 && IsPartCompatible(q->partition, partition, sumtParams.BitsLongsNeeded)==NO)
3372 q = &t->nodes[t->nNodes];
3373 q->partitionIndex = i;
3375 q->partition = partition;
3376 GetSummary(part->length, sumtParams.numRuns, part->count, &theStats, sumtParams.HPD);
3378 q->length = theStats.median;
3380 if (sumtParams.isClock == YES)
3382 GetSummary(part->height, sumtParams.numRuns, part->count, &theStats, sumtParams.HPD);
3383 q->depth = theStats.median;
3386 for (r = p->left; r!=NULL; r = r->sib)
3388 if (IsPartNested(r->partition, partition, sumtParams.BitsLongsNeeded) && r->depth >= q->depth)
3389 break; /* child is older then the node we try to add. Not good.*/
3391 if (p->depth <= q->depth)
3392 { /* New node older than the parent. Not good.*/
3393 r = p; /* Just to make r!=NULL*/
3397 if (sumtParams.isCalibrated == YES)
3399 GetSummary(part->age, sumtParams.numRuns, part->count, &theStats, sumtParams.HPD);
3400 q->age = theStats.median;
3403 for (r = p->left; r!=NULL; r = r->sib)
3405 if (freq < 1.00 && IsPartNested(r->partition, partition, sumtParams.BitsLongsNeeded) && r->age >= q->age)
3406 break; /* child is older then the node we try to add. Not good.*/
3408 if (p->age <= q->age)
3409 { /* New node older than the parent. Not good.*/
3410 r = p; /* Just to make r!=NULL*/
3415 if (r!=NULL && isFirstLoop)
3417 /* cancel the addition of the new node*/
3419 freqInterapted=freq;
3420 break; /* Finish creating the polytree */
3425 /* go through descendants of anc */
3427 for (r=p->left; r!=NULL; r=r ->sib)
3429 /* test if r is in the new partition or not */
3430 if ((r->x > 1 && IsPartNested(r->partition, partition, sumtParams.BitsLongsNeeded)) || (r->x == 1 && (partition[r->index / nBitsInALong] & (bitsLongOne << (r->index % nBitsInALong))) != 0))
3432 /* r is in the partition */
3442 /* r is not in the partition */
3450 /* terminate new sib-node chain */
3452 /* new node is last in old sib-node chain */
3458 /* singleton partition */
3460 if (nBits == sumtParams.numTaxa - 1)
3463 j = FirstTaxonInPartition(partition, sumtParams.BitsLongsNeeded); /* nbits == 1 */
3465 q->partitionIndex = i;
3466 q->partition = partition;
3467 numTerminalsEncountered++;
3468 GetSummary(part->length, sumtParams.numRuns, part->count, &theStats, sumtParams.HPD);
3469 q->length = theStats.median;
3470 if (sumtParams.isClock == YES)
3472 GetSummary(part->height, sumtParams.numRuns, part->count, &theStats, sumtParams.HPD);
3473 q->depth = theStats.median;
3474 if (q->anc->depth < q->depth)
3476 /* We never should get here because terminals always have 100% freq and they are younger than any other node that has 100% freq. */
3477 /* We should be careful with the trees with 0-brl generated under the fossilized birth-death prior! (<= is changed to <) */
3481 if (sumtParams.isCalibrated == YES)
3483 GetSummary(part->age, sumtParams.numRuns, part->count, &theStats, sumtParams.HPD);
3484 q->age = theStats.median;
3485 if (q->anc->age < q->age)
3487 /* We never should get here because terminals always have 100% freq and they are younger than any other node that has 100% freq. */
3488 /* We should be careful with the trees with 0-brl generated under the fossilized birth-death prior! (<= is changed to <) */
3501 goto treeConstruction;
3505 /* get downpass arrays */
3509 if (sumtParams.orderTaxa == YES)
3514 /* get downpass arrays */
3515 GetPolyDownPass(t2);
3518 if (sumtParams.orderTaxa == YES)
3522 /* draw tree to stdout and fp */
3523 MrBayesPrint ("\n%s Clade credibility values:\n\n", spacer);
3524 ShowConTree (stdout, t, 80, YES);
3525 if (logToFile == YES)
3526 ShowConTree (logFileFp, t, 80, YES);
3527 if (sumtParams.brlensDef == YES)
3529 MrBayesPrint ("\n");
3530 if (sumtParams.isClock == YES)
3531 MrBayesPrint ("%s Phylogram (based on median node depths):\n", spacer);
3533 MrBayesPrint ("%s Phylogram (based on average branch lengths):\n", spacer);
3536 MrBayesPrint ("%s Warning. Phylogram containing all nodes with credibility values exceeding\n",spacer);
3537 MrBayesPrint ("%s the level set by Contype could not be constructed.\n",spacer);
3538 MrBayesPrint ("%s Only nodes with credibility values exceeding %.2f%% (percentage of trees\n", spacer, freqInterapted*100);
3539 MrBayesPrint ("%s where the node is present) were included in the phylogram.\n", spacer);
3541 MrBayesPrint ("\n");
3542 ShowConPhylogram (stdout, t2, 80);
3543 if (logToFile == YES)
3544 ShowConPhylogram (logFileFp, t2, 80);
3547 /* print taxa block */
3548 MrBayesPrintf (fpCon, "begin taxa;\n");
3549 MrBayesPrintf (fpCon, "\tdimensions ntax=%d;\n", sumtParams.numTaxa);
3550 MrBayesPrintf (fpCon, "\ttaxlabels\n", sumtParams.numTaxa);
3551 for (i=0; i<sumtParams.numTaxa; i++)
3553 for (j=0; j<t2->nNodes; j++)
3554 if (t2->nodes[j].index == i)
3556 MrBayesPrintf (fpCon, "\t\t%s\n", t2->nodes[j].label);
3558 MrBayesPrintf (fpCon, "\t\t;\nend;\n");
3560 MrBayesPrintf (fpCon, "begin trees;\n");
3561 MrBayesPrintf (fpCon, "\ttranslate\n");
3562 for (i=0; i<sumtParams.numTaxa; i++)
3564 for (j=0; j<t2->nNodes; j++)
3565 if (t2->nodes[j].index == i)
3567 if (i == sumtParams.numTaxa-1)
3568 MrBayesPrintf (fpCon, "\t\t%d\t%s\n", t2->nodes[i].index+1, t2->nodes[i].label);
3570 MrBayesPrintf (fpCon, "\t\t%d\t%s,\n", t2->nodes[i].index+1, t2->nodes[i].label);
3572 MrBayesPrintf (fpCon, "\t\t;\n");
3573 if (sumtParams.consensusFormat == SIMPLE)
3574 PrintConTree(fpCon, t2);
3575 else if (sumtParams.consensusFormat == FIGTREE)
3576 PrintFigTreeConTree(fpCon, t2, treeParts);
3577 MrBayesPrintf (fpCon, "end;\n");
3587 MrBayesPrint ("%d", numTreeParts); /* just because I am tired of seeing the unused parameter error msg */
3591 MrBFlt CppEvolRate (PolyTree *t, PolyNode *p, int eSet)
3594 MrBFlt ancRate, branchRate, *rate, *pos;
3597 nEvents = t->nEvents[eSet][p->index];
3598 pos = t->position[eSet][p->index];
3599 rate = t->rateMult[eSet][p->index];
3601 /* note that event positions are from top of branch (more recent, descendant tip) */
3603 // if (t->eType[eSet] == CPPm)
3605 for (q=p; q->anc != NULL; q=q->anc)
3607 for (i=0; i<t->nEvents[eSet][p->index]; i++)
3608 ancRate *= t->rateMult[eSet][p->index][i];
3612 branchRate = rate[0] * pos[0];
3613 for (i=1; i<nEvents; i++)
3615 branchRate += (pos[i] - pos[i-1]);
3616 branchRate *= rate[i];
3618 branchRate += 1.0 - pos[nEvents-1];
3619 branchRate *= ancRate;
3622 branchRate = ancRate;
3625 else if (t->eType[eSet] == CPPi)
3627 for (q=p; q->anc != NULL; q=q->anc)
3629 if (t->nEvents[eSet][p->index]>0)
3631 ancRate = t->rateMult[eSet][p->index][0];
3637 branchRate = ancRate * (1.0 - pos[nEvents-1]);
3638 for (i=nEvents-2; i>=0; i--)
3640 branchRate += (rate[i+1] * (pos[i+1] - pos[i]));
3642 branchRate += (rate[0] * pos[0]);
3645 branchRate = ancRate;
3653 int DoCompareTree (void)
3655 int i, j, k, n, longestLineLength, brlensDef[2], numTreesInLastBlock[2],
3656 lastTreeBlockBegin[2], lastTreeBlockEnd[2], xaxis, yaxis, starHolder[80],
3657 minNumTrees, screenWidth, screenHeigth, numY[60], nSamples;
3658 RandLong temporarySeed;
3661 MrBFlt xProb, yProb, xInc, yInc, xUpper, xLower, yUpper, yLower, *dT1=NULL, *dT2=NULL, *dT3=NULL, d1, d2, d3,
3662 meanY[60], xVal, yVal, minX, minY, maxX, maxY, sums[3];
3663 char *s=NULL, prCh, treeName[2][100];
3666 PartCtr **treeParts=NULL;
3667 Tree *tree1=NULL, *tree2=NULL;
3668 SumtFileInfo sumtFileInfo;
3670 # if defined (MPI_ENABLED)
3675 /* Make sure we read trees using DoSumtTree() code instead of with the user tree code */
3676 inComparetreeCommand = YES;
3678 /* set file pointer to NULL */
3681 strcpy(treeName[0],"tree"); //in case if parameter is not specified in a .t file
3682 strcpy(treeName[1],"tree");
3684 /* Check that a data set has been read in. We check taxon names against
3686 if (isTaxsetDef == NO)
3688 MrBayesPrint ("%s A matrix or set of taxon labels must be specified before comparetree can be used\n", spacer);
3692 /* open output files for summary information (two files); check if we want to overwrite previous results */
3693 if (OpenComptFiles () == ERROR)
3696 MrBayesPrint ("%s Examining files ...\n", spacer);
3698 /* Examine first file */
3699 if (ExamineSumtFile(comptreeParams.comptFileName1, &sumtFileInfo, treeName[0], &(brlensDef[0])) == ERROR)
3703 longestLineLength = sumtFileInfo.longestLineLength;
3704 numTreesInLastBlock[0] = sumtFileInfo.numTreesInLastBlock;
3705 lastTreeBlockBegin[0] = sumtFileInfo.lastTreeBlockBegin;
3706 lastTreeBlockEnd[0] = sumtFileInfo.lastTreeBlockEnd;
3708 /* Examine second file */
3709 if (ExamineSumtFile(comptreeParams.comptFileName2, &sumtFileInfo, treeName[1], &brlensDef[1]) == ERROR)
3713 if (longestLineLength < sumtFileInfo.longestLineLength)
3714 longestLineLength = sumtFileInfo.longestLineLength;
3715 numTreesInLastBlock[1] = sumtFileInfo.numTreesInLastBlock;
3716 lastTreeBlockBegin[1] = sumtFileInfo.lastTreeBlockBegin;
3717 lastTreeBlockEnd[1] = sumtFileInfo.lastTreeBlockEnd;
3719 /* Check whether we should work with brlens */
3720 if (brlensDef[0] == YES && brlensDef[1] == YES)
3721 sumtParams.brlensDef = YES;
3723 sumtParams.brlensDef = NO;
3725 /* Allocate space for command string */
3726 longestLineLength += 10;
3727 s = (char *)SafeMalloc((size_t)longestLineLength * sizeof(char));
3730 MrBayesPrint ("%s Problem allocating string for reading tree file\n", spacer);
3734 /* Allocate space for packed trees */
3735 if (chainParams.relativeBurnin == YES)
3737 numPackedTrees[0] = numTreesInLastBlock[0] - (int)(chainParams.burninFraction * numTreesInLastBlock[0]);
3738 numPackedTrees[1] = numTreesInLastBlock[1] - (int)(chainParams.burninFraction * numTreesInLastBlock[1]);
3742 numPackedTrees[0] = numTreesInLastBlock[0] - chainParams.chainBurnIn;
3743 numPackedTrees[1] = numTreesInLastBlock[1] - chainParams.chainBurnIn;
3745 if (memAllocs[ALLOC_PACKEDTREES] == YES)
3747 MrBayesPrint ("%s packedTreeList is already allocated\n", spacer);
3750 packedTreeList[0] = (PackedTree *) SafeCalloc(numPackedTrees[0]+numPackedTrees[1], sizeof(PackedTree));
3751 packedTreeList[1] = packedTreeList[0] + numPackedTrees[0];
3752 if (!packedTreeList[0])
3754 MrBayesPrint ("%s Problem allocating packed tree list\n", spacer);
3757 memAllocs[ALLOC_PACKEDTREES] = YES;
3759 /* Tell user we are ready to go */
3760 MrBayesPrint ("%s Summarizing trees in files \"%s\" and \"%s\"\n", spacer,
3761 comptreeParams.comptFileName1,
3762 comptreeParams.comptFileName2);
3764 if (chainParams.relativeBurnin == YES)
3765 MrBayesPrint ("%s Using relative burnin ('relburnin=yes'), discarding the first %.0f %% ('burninfrac=%1.2f') of sampled trees\n",
3766 spacer, chainParams.burninFraction*100.0, chainParams.burninFraction);
3768 MrBayesPrint ("%s Using absolute burnin ('relburnin=no'), discarding the first %d ('burnin=%d') sampled trees\n",
3769 spacer, chainParams.chainBurnIn, chainParams.chainBurnIn);
3771 MrBayesPrint ("%s Writing statistics to file %s.<dists|pairs>\n", spacer, comptreeParams.comptOutfile);
3773 /* Set up cheap status bar. */
3774 MrBayesPrint ("\n%s Tree reading status:\n\n", spacer);
3775 MrBayesPrint ("%s 0 10 20 30 40 50 60 70 80 90 100\n", spacer);
3776 MrBayesPrint ("%s v-------v-------v-------v-------v-------v-------v-------v-------v-------v-------v\n", spacer);
3777 MrBayesPrint ("%s *", spacer);
3780 /* Read file 1 for real */
3781 if ((fp = OpenTextFileR(comptreeParams.comptFileName1)) == NULL)
3784 /* ...and fast forward to beginning of last tree block (skipping begin trees). */
3785 for (i=0; i<lastTreeBlockBegin[0] + 1; i++)
3787 if (fgets (s, longestLineLength, fp) == NULL)
3789 printf ("Error in function: %s at line: %d in file: %s", __FUNCTION__, __LINE__, __FILE__);
3793 /* Calculate burnin */
3794 if (chainParams.relativeBurnin == YES)
3795 comptreeParams.burnin = (int)(chainParams.burninFraction * numTreesInLastBlock[0]);
3797 comptreeParams.burnin = chainParams.chainBurnIn;
3799 /* Initialize sumtParams struct */
3800 numUniqueSplitsFound = numUniqueTreesFound = 0;
3801 sumtParams.runId = 0;
3802 strcpy(sumtParams.curFileName, comptreeParams.comptFileName1);
3803 sumtParams.tree = AllocatePolyTree (numTaxa);
3804 AllocatePolyTreePartitions (sumtParams.tree);
3805 sumtParams.numTreesEncountered = sumtParams.numTreesSampled = 0;
3806 sumtParams.numFileTrees = (int *) SafeCalloc (2*2+2*numTaxa, sizeof(int));
3807 sumtParams.numFileTreesSampled = sumtParams.numFileTrees + sumtParams.numRuns;
3808 sumtParams.order = sumtParams.numFileTrees + 2*sumtParams.numRuns;
3809 sumtParams.absentTaxa = sumtParams.numFileTrees + 2*sumtParams.numRuns + numTaxa;
3810 sumtParams.numTreesInLastBlock = numTreesInLastBlock[0];
3811 if (!sumtParams.numFileTrees)
3813 MrBayesPrint ("%s Problems allocating sumtParams.numFileTrees in DoSumt()\n", spacer);
3817 memAllocs[ALLOC_SUMTPARAMS] = YES;
3819 /* ... and parse the file */
3820 expecting = Expecting(COMMAND);
3822 ResetTranslateTable();
3823 for (i=0; i<lastTreeBlockEnd[0] - lastTreeBlockBegin[0] - 1; i++)
3825 if (fgets (s, longestLineLength, fp) == NULL)
3827 printf ("Error in function: %s at line: %d in file: %s", __FUNCTION__, __LINE__, __FILE__);
3829 /*MrBayesPrint ("%s", s);*/
3830 if (ParseCommand (s) == ERROR)
3834 ResetTranslateTable();
3836 /* Check that at least one tree was read in. */
3837 if (sumtParams.numFileTreesSampled[0] <= 0)
3839 MrBayesPrint ("%s No trees read in\n", spacer);
3843 /* ... and close file */
3846 /* Read file 2 for real */
3847 if ((fp = OpenTextFileR(comptreeParams.comptFileName2)) == NULL)
3850 /* ...and fast forward to beginning of last tree block. */
3851 for (i=0; i<lastTreeBlockBegin[1] + 1; i++)
3853 if (fgets (s, longestLineLength, fp) == NULL)
3855 printf ("Error in function: %s at line: %d in file: %s", __FUNCTION__, __LINE__, __FILE__);
3859 /* Renitialize sumtParams struct */
3860 sumtParams.runId = 1;
3861 strcpy (sumtParams.curFileName, comptreeParams.comptFileName2);
3863 /* Calculate burnin */
3864 if (chainParams.relativeBurnin == YES)
3865 comptreeParams.burnin = (int)(chainParams.burninFraction * numTreesInLastBlock[1]);
3867 comptreeParams.burnin = chainParams.chainBurnIn;
3869 /* ... and parse the file */
3870 expecting = Expecting(COMMAND);
3872 ResetTranslateTable();
3873 for (i=0; i<lastTreeBlockEnd[1] - lastTreeBlockBegin[1] - 1; i++)
3875 if (fgets (s, longestLineLength, fp) == NULL)
3877 printf ("Error in function: %s at line: %d in file: %s", __FUNCTION__, __LINE__, __FILE__);
3879 /*MrBayesPrint ("%s", s);*/
3880 if (ParseCommand (s) == ERROR)
3884 ResetTranslateTable();
3886 /* Check that at least one tree was read in. */
3887 if (sumtParams.numFileTreesSampled[1] <= 0)
3889 MrBayesPrint ("%s No trees read in\n", spacer);
3893 /* ... and close file */
3896 /* Now finish cheap status bar. */
3897 if (numAsterices < 80)
3898 for (i=0; i<80 - numAsterices; i++)
3900 MrBayesPrint ("\n\n");
3902 /* tell user how many trees were successfully read */
3903 MrBayesPrint ("%s Read %d trees from last tree block of file \"%s\" (sampling %d of them)\n", spacer,
3904 sumtParams.numFileTrees[0],
3905 comptreeParams.comptFileName1,
3906 sumtParams.numFileTreesSampled[0]);
3907 MrBayesPrint ("%s Read %d trees from last tree block of file \"%s\" (sampling %d of them)\n", spacer,
3908 sumtParams.numFileTrees[1],
3909 comptreeParams.comptFileName2,
3910 sumtParams.numFileTreesSampled[1]);
3912 /* Extract partition counter pointers */
3913 treeParts = (PartCtr **) SafeCalloc ((size_t)(numUniqueSplitsFound), sizeof(PartCtr *));
3915 PartCtrUppass(partCtrRoot, treeParts, &i);
3917 /* Sort taxon partitions (clades, splits) ... */
3918 SortPartCtr (treeParts, 0, numUniqueSplitsFound-1);
3920 /* print to screen */
3921 MrBayesPrint (" \n");
3922 MrBayesPrint ("%s General explanation: \n", spacer);
3923 MrBayesPrint (" \n");
3924 MrBayesPrint ("%s In an unrooted tree, a taxon bipartition (split) is specified by removing a \n", spacer);
3925 MrBayesPrint ("%s branch, thereby dividing the species into those to the left and those to the \n", spacer);
3926 MrBayesPrint ("%s right of the branch. Here, taxa to one side of the removed branch are denoted \n", spacer);
3927 MrBayesPrint ("%s '.' and those to the other side are denoted '*'. Specifically, the '.' symbol \n", spacer);
3928 MrBayesPrint ("%s is used for the taxa on the same side as the outgroup. \n", spacer);
3929 MrBayesPrint (" \n");
3930 MrBayesPrint ("%s In a rooted or clock tree, the '*' symbol is simply used to denote the taxa \n", spacer);
3931 MrBayesPrint ("%s that are included in a particular group (clade), that is, all the descendants \n", spacer);
3932 MrBayesPrint ("%s of a particular branch in the tree. Taxa that are not included are denoted \n", spacer);
3933 MrBayesPrint ("%s using the '.' symbol. \n", spacer);
3934 MrBayesPrint (" \n");
3935 MrBayesPrint ("%s The output includes the ID of the encountered clades or splits (sorted from \n", spacer);
3936 MrBayesPrint ("%s highest to lowest probability), the bipartition or clade in '.*' format, \n", spacer);
3937 MrBayesPrint ("%s number of times the bipartition or clade was observed in the first tree file \n", spacer);
3938 MrBayesPrint ("%s served in the first tree file, the number of times the bipartition was, \n", spacer);
3939 MrBayesPrint ("%s observed in the second tree file, the proportion of the time the bipartition \n", spacer);
3940 MrBayesPrint ("%s was found in the first tree file, and the proportion of the time the bi- \n", spacer);
3941 MrBayesPrint ("%s partition was found in the second tree file. \n", spacer);
3942 MrBayesPrint (" \n");
3943 MrBayesPrint ("%s List of taxa in bipartitions: \n", spacer);
3944 MrBayesPrint (" \n");
3946 for (k=0; k<numTaxa; k++)
3948 if (sumtParams.absentTaxa[k] == NO && taxaInfo[k].isDeleted == NO)
3950 MrBayesPrint ("%s %4d -- %s\n", spacer, j++, taxaNames[k]);
3953 MrBayesPrint (" \n");
3954 MrBayesPrint ("%s List of taxon bipartitions found in tree files: \n\n", spacer);
3956 i = (int)(log10(sumtParams.numTreesSampled)) - 1;
3959 j = sumtParams.numTaxa - 8;
3962 MrBayesPrint ("%s ID -- Partition%*c No1%*c No2%*c Freq1 Freq2\n",
3963 spacer, j, ' ', i, ' ', i, ' ');
3965 mask = SafeCalloc (sumtParams.BitsLongsNeeded, sizeof(BitsLong));
3966 for (i=0; i<sumtParams.numTaxa; i++)
3968 for (i=0; i<numUniqueSplitsFound; i++)
3971 if (IsBitSet(localOutGroup, x->partition) == YES && sumtParams.isRooted == NO)
3972 FlipBits(x->partition, sumtParams.BitsLongsNeeded, mask);
3973 if ((MrBFlt)x->totCount/(MrBFlt)sumtParams.numTreesSampled >= comptreeParams.minPartFreq)
3975 MrBayesPrint ("%s %4d -- ", spacer, i+1);
3976 ShowParts (stdout, x->partition, sumtParams.numTaxa);
3978 j = (int)(log10(sumtParams.numTreesSampled)) + 1;
3981 MrBayesPrint (" %*d %*d %1.3lf %1.3lf\n",
3982 j, x->count[0], j, x->count[1],
3983 (MrBFlt)x->count[0]/(MrBFlt)sumtParams.numFileTreesSampled[0],
3984 (MrBFlt)x->count[1]/(MrBFlt)sumtParams.numFileTreesSampled[1]);
3986 MrBayesPrintf (fpParts, "%d\t%d\t%d\t%1.3lf\t%1.3lf\n",
3987 i+1, x->count[0], x->count[1],
3988 (MrBFlt)x->count[0]/(MrBFlt)sumtParams.numFileTreesSampled[0],
3989 (MrBFlt)x->count[1]/(MrBFlt)sumtParams.numFileTreesSampled[1]);
3994 /* make a nifty graph plotting frequencies of clades found in the two tree files */
3995 MrBayesPrint (" \n");
3996 MrBayesPrint ("%s Bivariate plot of clade probabilities: \n", spacer);
3997 MrBayesPrint (" \n");
3998 MrBayesPrint ("%s This graph plots the probabilities of clades found in file 1 (the x-axis) \n", spacer);
3999 MrBayesPrint ("%s against the probabilities of the same clades found in file 2 (the y-axis). \n", spacer);
4000 MrBayesPrint (" \n");
4001 xInc = (MrBFlt) (1.0 / 80.0);
4002 yInc = (MrBFlt) (1.0 / 40.0);
4004 yLower = yUpper - yInc;
4005 for (yaxis=39; yaxis>=0; yaxis--)
4008 xUpper = xLower + xInc;
4009 for (xaxis=0; xaxis<80; xaxis++)
4011 starHolder[xaxis] = 0;
4012 for (i=0; i<numUniqueSplitsFound; i++)
4015 xProb = (MrBFlt)x->count[0]/(MrBFlt)sumtParams.numFileTreesSampled[0];
4016 yProb = (MrBFlt)x->count[1]/(MrBFlt)sumtParams.numFileTreesSampled[1];
4017 if ((xProb > xLower || (xProb == 0.0 && xaxis == 0)) && (xProb <= xUpper || (xProb == 1.0 && xaxis == 79))
4018 && (yProb > yLower || (yProb == 0.0 && yaxis == 0)) && (yProb <= yUpper || (yProb == 1.0 && yaxis == 39)))
4019 starHolder[xaxis] = 1;
4022 xUpper = xLower + xInc;
4025 MrBayesPrint ("%s ", spacer);
4026 for (xaxis=0; xaxis<80; xaxis++)
4029 if ((xaxis == 0 && yaxis == 0) || (xaxis == 79 && yaxis == 39))
4031 else if ((xaxis == 0 && yaxis == 39) || (xaxis == 79 && yaxis == 0))
4033 else if ((yaxis == 0 || yaxis == 39) && xaxis > 0 && xaxis < 79)
4035 else if ((xaxis == 0 || xaxis == 79) && yaxis > 0 && yaxis < 39)
4037 if (starHolder[xaxis] == 1)
4039 MrBayesPrint ("%c", prCh);
4042 MrBayesPrint (" 1.00\n");
4043 else if (yaxis == 0)
4044 MrBayesPrint (" 0.00\n");
4046 MrBayesPrint ("\n");
4049 yLower = yUpper - yInc;
4052 MrBayesPrint ("%s ^ ^\n", spacer);
4053 MrBayesPrint ("%s 0.00 1.00\n", spacer);
4055 /* get tree-to-tree distances: first allocate some space */
4056 minNumTrees = sumtParams.numFileTreesSampled[0];
4057 if (sumtParams.numFileTreesSampled[1] < minNumTrees)
4058 minNumTrees = sumtParams.numFileTreesSampled[1];
4059 dT1 = (MrBFlt *) SafeMalloc (3 * (size_t)minNumTrees * sizeof(MrBFlt));
4060 tree1 = AllocateFixedTree (sumtParams.numTaxa, sumtParams.isRooted);
4061 tree2 = AllocateFixedTree (sumtParams.numTaxa, sumtParams.isRooted);
4062 if (!dT1 || !tree1 || !tree2)
4064 MrBayesPrint ("%s Problem allocating topological distances\n", spacer);
4067 dT2 = dT1 + minNumTrees;
4068 dT3 = dT2 + minNumTrees;
4070 for (i=0; i<minNumTrees; i++)
4072 if (sumtParams.isRooted == NO)
4074 RetrieveUTree (tree1, packedTreeList[0][i].order, packedTreeList[0][i].brlens);
4075 RetrieveUTree (tree2, packedTreeList[1][i].order, packedTreeList[1][i].brlens);
4079 RetrieveRTree (tree1, packedTreeList[0][i].order, packedTreeList[0][i].brlens);
4080 RetrieveRTree (tree2, packedTreeList[1][i].order, packedTreeList[1][i].brlens);
4082 /* Allocate and set partitions now that we have a tree */
4085 AllocateTreePartitions(tree1);
4086 AllocateTreePartitions(tree2);
4090 ResetTreePartitions(tree1);
4091 ResetTreePartitions(tree2);
4093 CalculateTreeToTreeDistance (tree1, tree2, &d1, &d2, &d3);
4099 for (i=0; i<minNumTrees; i++)
4101 /*MrBayesPrint ("%s %4d -- %lf %lf %lf\n", spacer, i+1, dT1[i], dT2[i], dT3[i]);*/
4102 if (sumtParams.brlensDef == YES)
4103 MrBayesPrintf (fpDists, "%d\t%lf\t%lf\t%lf\n", i+1, dT1[i], dT2[i], dT3[i]);
4105 MrBayesPrintf (fpDists, "%d\t%lf\n", i+1, dT1[i]);
4108 /* print x-y plot of log likelihood vs. generation */
4109 MrBayesPrint ("\n");
4110 MrBayesPrint ("%s Rough plots of generation (x-axis) versus the measures of tree- \n", spacer);
4111 MrBayesPrint ("%s to-tree distances (y-axis). \n", spacer);
4112 MrBayesPrint ("\n");
4113 MrBayesPrint ("%s Distance(Robinson-Foulds):\n", spacer);
4114 MrBayesPrint ("\n");
4115 screenWidth = 60; /* don't change this without changing numY and meanY, declared above */
4117 minX = minY = 1000000000.0;
4118 maxX = maxY = -1000000000.0;
4119 for (i=0; i<minNumTrees; i++)
4121 xVal = (MrBFlt) (i + comptreeParams.burnin);
4132 for (i=0; i<screenWidth; i++)
4137 for (i=0; i<minNumTrees; i++)
4139 xVal = (MrBFlt) (i + comptreeParams.burnin);
4141 k = (int)(((xVal - minX) / (maxX - minX)) * screenWidth);
4142 if (k >= screenWidth)
4143 k = screenWidth - 1;
4147 MrBayesPrint ("\n%s +", spacer);
4148 for (i=0; i<screenWidth; i++)
4150 MrBayesPrint ("+ %1.2lf\n", maxY);
4151 for (j=screenHeigth-1; j>=0; j--)
4153 MrBayesPrint ("%s |", spacer);
4154 for (i=0; i<screenWidth; i++)
4158 if (meanY[i] / numY[i] > (((maxY - minY)/screenHeigth)*j)+minY && meanY[i] / numY[i] <= (((maxY - minY)/screenHeigth)*(j+1))+minY)
4168 MrBayesPrint ("|\n");
4170 MrBayesPrint ("%s +", spacer);
4171 for (i=0; i<screenWidth; i++)
4173 if (numY[i] > 0 && meanY[i] / numY[i] <= minY)
4175 else if (i % (screenWidth/10) == 0 && i != 0)
4180 MrBayesPrint ("+ %1.2lf\n", minY);
4181 MrBayesPrint ("%s ^", spacer);
4182 for (i=0; i<screenWidth; i++)
4184 MrBayesPrint ("^\n");
4185 MrBayesPrint ("%s %1.0lf", spacer, minX);
4186 for (i=0; i<screenWidth; i++)
4188 MrBayesPrint ("%1.0lf\n\n", maxX);
4190 if (sumtParams.brlensDef == YES)
4194 MrBayesPrint ("\n");
4196 MrBayesPrint ("%s Distance(Robinson-Foulds with branch lengths):\n", spacer);
4198 MrBayesPrint ("%s Distance(Robinson-Foulds with scaled branch lengths):\n", spacer);
4199 MrBayesPrint ("\n");
4200 screenWidth = 60; /* don't change this without changing numY and meanY, declared above */
4202 minX = minY = 1000000000.0;
4203 maxX = maxY = -1000000000.0;
4204 for (i=0; i<minNumTrees; i++)
4206 xVal = (MrBFlt) (i + comptreeParams.burnin);
4220 for (i=0; i<screenWidth; i++)
4225 for (i=0; i<minNumTrees; i++)
4227 xVal = (MrBFlt) (i + comptreeParams.burnin);
4232 k = (int)(((xVal - minX) / (maxX - minX)) * screenWidth);
4233 if (k >= screenWidth)
4234 k = screenWidth - 1;
4238 MrBayesPrint ("\n%s +", spacer);
4239 for (i=0; i<screenWidth; i++)
4241 MrBayesPrint ("+ %1.2lf\n", maxY);
4242 for (j=screenHeigth-1; j>=0; j--)
4244 MrBayesPrint ("%s |", spacer);
4245 for (i=0; i<screenWidth; i++)
4249 if (meanY[i] / numY[i] > (((maxY - minY)/screenHeigth)*j)+minY && meanY[i] / numY[i] <= (((maxY - minY)/screenHeigth)*(j+1))+minY)
4259 MrBayesPrint ("|\n");
4261 MrBayesPrint ("%s +", spacer);
4262 for (i=0; i<screenWidth; i++)
4264 if (numY[i] > 0 && meanY[i] / numY[i] <= minY)
4266 else if (i % (screenWidth/10) == 0 && i != 0)
4271 MrBayesPrint ("+ %1.2lf\n", minY);
4272 MrBayesPrint ("%s ^", spacer);
4273 for (i=0; i<screenWidth; i++)
4275 MrBayesPrint ("^\n");
4276 MrBayesPrint ("%s %1.0lf", spacer, minX);
4277 for (i=0; i<screenWidth; i++)
4279 MrBayesPrint ("%1.0lf\n\n", maxX);
4283 /* calculate average tree-to-tree distances */
4284 curTime = time(NULL);
4285 temporarySeed = (RandLong)curTime;
4286 if (temporarySeed < 0)
4287 temporarySeed = -temporarySeed;
4288 sums[0] = sums[1] = sums[2] = 0.0;
4290 for (n=0; n<nSamples; n++)
4292 i = (int) RandomNumber(&temporarySeed) * minNumTrees;
4293 j = (int) RandomNumber(&temporarySeed) * minNumTrees;
4294 if (sumtParams.isRooted == NO)
4296 RetrieveUTree (tree1, packedTreeList[0][i].order, packedTreeList[0][i].brlens);
4297 RetrieveUTree (tree2, packedTreeList[1][j].order, packedTreeList[1][j].brlens);
4299 else /* if (sumtParams.isRooted == YES) */
4301 RetrieveRTree (tree1, packedTreeList[0][i].order, packedTreeList[0][i].brlens);
4302 RetrieveRTree (tree2, packedTreeList[1][j].order, packedTreeList[1][j].brlens);
4304 CalculateTreeToTreeDistance (tree1, tree2, &d1, &d2, &d3);
4309 MrBayesPrint ("%s Mean tree-to-tree distances, based on %d trees randomly sampled from both files:\n\n", spacer, nSamples);
4310 MrBayesPrint ("%s Mean(Robinson-Foulds) = %1.3lf\n", spacer, sums[0]/nSamples);
4311 if (sumtParams.brlensDef == YES)
4313 MrBayesPrint ("%s Mean(Robinson-Foulds with branch lengths) = %1.3lf\n", spacer, sums[1]/nSamples);
4314 MrBayesPrint ("%s Mean(Robinson-Foulds with scaled branch lengths) = %1.3lf\n", spacer, sums[2]/nSamples);
4316 MrBayesPrint ("\n");
4318 /* free memory and file pointers */
4323 if (memAllocs[ALLOC_PACKEDTREES] == YES)
4325 for (i=0; i<numPackedTrees[0]+numPackedTrees[1]; i++)
4327 free(packedTreeList[0][i].order);
4328 free(packedTreeList[0][i].brlens);
4330 free(packedTreeList[0]);
4331 memAllocs[ALLOC_PACKEDTREES] = NO;
4334 /* free sumtParams */
4339 SafeFclose (&fpParts);
4340 SafeFclose (&fpDists);
4342 /* free pointer array to partitions, part and tree counters */
4344 FreePartCtr (partCtrRoot);
4345 FreeTreeCtr (treeCtrRoot);
4349 /* reset taxon set */
4352 # if defined (MPI_ENABLED)
4356 expecting = Expecting(COMMAND);
4357 inComparetreeCommand = NO;
4364 /* free sumtParams */
4370 if (memAllocs[ALLOC_PACKEDTREES] == YES)
4372 for (i=0; i<numPackedTrees[0]+numPackedTrees[1]; i++)
4374 free(packedTreeList[0][i].order);
4375 free(packedTreeList[0][i].brlens);
4377 free(packedTreeList[0]);
4378 memAllocs[ALLOC_PACKEDTREES] = NO;
4381 /* reset taxon set */
4384 /* free pointer array to partitions, part and tree counters */
4386 FreePartCtr (partCtrRoot);
4387 FreeTreeCtr (treeCtrRoot);
4392 SafeFclose (&fpParts);
4393 SafeFclose (&fpDists);
4394 strcpy (spacer, "");
4396 expecting = Expecting(COMMAND);
4397 inComparetreeCommand = NO;
4403 int DoCompareTreeParm (char *parmName, char *tkn)
4409 if (expecting == Expecting(PARAMETER))
4411 expecting = Expecting(EQUALSIGN);
4415 if (!strcmp(parmName, "Xxxxxxxxxx"))
4417 expecting = Expecting(PARAMETER);
4418 expecting |= Expecting(SEMICOLON);
4420 /* set Filename (comptreeParams.comptFileName1) *************************************************/
4421 else if (!strcmp(parmName, "Filename1"))
4423 if (expecting == Expecting(EQUALSIGN))
4425 expecting = Expecting(ALPHA);
4428 else if (expecting == Expecting(ALPHA))
4430 strcpy (comptreeParams.comptFileName1, tkn);
4431 MrBayesPrint ("%s Setting comparetree filename 1 to %s\n", spacer, comptreeParams.comptFileName1);
4432 expecting = Expecting(PARAMETER) | Expecting(SEMICOLON);
4437 /* set Filename (comptreeParams.comptFileName2) *************************************************/
4438 else if (!strcmp(parmName, "Filename2"))
4440 if (expecting == Expecting(EQUALSIGN))
4442 expecting = Expecting(ALPHA);
4445 else if (expecting == Expecting(ALPHA))
4447 strcpy (comptreeParams.comptFileName2, tkn);
4448 MrBayesPrint ("%s Setting comparetree filename 2 to %s\n", spacer, comptreeParams.comptFileName2);
4449 expecting = Expecting(PARAMETER) | Expecting(SEMICOLON);
4454 /* set Filename (comptreeParams.comptOutfile) ***************************************************/
4455 else if (!strcmp(parmName, "Outputname"))
4457 if (expecting == Expecting(EQUALSIGN))
4459 expecting = Expecting(ALPHA);
4462 else if (expecting == Expecting(ALPHA))
4464 strcpy (comptreeParams.comptOutfile, tkn);
4465 MrBayesPrint ("%s Setting comparetree output file to %s\n", spacer, comptreeParams.comptOutfile);
4466 expecting = Expecting(PARAMETER) | Expecting(SEMICOLON);
4471 /* set Relburnin (chainParams.relativeBurnin) ********************************************************/
4472 else if (!strcmp(parmName, "Relburnin"))
4474 if (expecting == Expecting(EQUALSIGN))
4475 expecting = Expecting(ALPHA);
4476 else if (expecting == Expecting(ALPHA))
4478 if (IsArgValid(tkn, tempStr) == NO_ERROR)
4480 if (!strcmp(tempStr, "Yes"))
4481 chainParams.relativeBurnin = YES;
4483 chainParams.relativeBurnin = NO;
4487 MrBayesPrint ("%s Invalid argument for Relburnin\n", spacer);
4490 if (chainParams.relativeBurnin == YES)
4491 MrBayesPrint ("%s Using relative burnin (a fraction of samples discarded).\n", spacer);
4493 MrBayesPrint ("%s Using absolute burnin (a fixed number of samples discarded).\n", spacer);
4494 expecting = Expecting(PARAMETER) | Expecting(SEMICOLON);
4502 /* set Burnin (chainParams.chainBurnIn) *******************************************************/
4503 else if (!strcmp(parmName, "Burnin"))
4505 if (expecting == Expecting(EQUALSIGN))
4506 expecting = Expecting(NUMBER);
4507 else if (expecting == Expecting(NUMBER))
4509 sscanf (tkn, "%d", &tempI);
4510 chainParams.chainBurnIn = tempI;
4511 MrBayesPrint ("%s Setting burnin to %d\n", spacer, chainParams.chainBurnIn);
4512 expecting = Expecting(PARAMETER) | Expecting(SEMICOLON);
4517 /* set Burninfrac (chainParams.burninFraction) ************************************************************/
4518 else if (!strcmp(parmName, "Burninfrac"))
4520 if (expecting == Expecting(EQUALSIGN))
4521 expecting = Expecting(NUMBER);
4522 else if (expecting == Expecting(NUMBER))
4524 sscanf (tkn, "%lf", &tempD);
4527 MrBayesPrint ("%s Burnin fraction too low (< 0.01)\n", spacer);
4532 MrBayesPrint ("%s Burnin fraction too high (> 0.50)\n", spacer);
4535 chainParams.burninFraction = tempD;
4536 MrBayesPrint ("%s Setting burnin fraction to %.2f\n", spacer, chainParams.burninFraction);
4537 expecting = Expecting(PARAMETER) | Expecting(SEMICOLON);
4544 /* set Minpartfreq (comptreeParams.minPartFreq) *******************************************************/
4545 else if (!strcmp(parmName, "Minpartfreq"))
4547 if (expecting == Expecting(EQUALSIGN))
4548 expecting = Expecting(NUMBER);
4549 else if (expecting == Expecting(NUMBER))
4551 sscanf (tkn, "%lf", &tempD);
4552 comptreeParams.minPartFreq = tempD;
4553 MrBayesPrint ("%s Including partitions with probability greater than or equal to %lf in summary statistics\n", spacer, comptreeParams.minPartFreq);
4554 expecting = Expecting(PARAMETER) | Expecting(SEMICOLON);
4567 int DoCompRefTree (void)
4569 /* Compare a tree file with the reference tree files to generate the SDSFs.
4570 Use parameters in CompareTree and MCMCP (lazy option) */
4572 char outName[130], inName[130], inRefName[130], treeName[100], *lineBuf=NULL, *s;
4573 FILE *fpTre=NULL, *fpOut=NULL;
4574 int i, n, longestL=0, burnin, gen, nRefRun, nTre[100]={0};
4575 SumtFileInfo tFileInfo;
4578 # if defined (MPI_ENABLED)
4583 /* Check that a data set has been read in. We check taxon names against those read in. */
4584 if (isTaxsetDef == NO)
4586 MrBayesPrint ("%s A matrix or set of taxon labels must be specified before compareref can be used\n", spacer);
4590 /* this is a hack to account for the additional comparing tree sample
4591 chainParams.numRuns is used in AddTreeToPartitionCounters to alloc mem correctly */
4592 nRefRun = chainParams.numRuns;
4593 chainParams.numRuns += 1;
4596 if ((t = AllocateTree (numLocalTaxa)) == NULL)
4598 MrBayesPrint ("%s Problem allocating diagn tree\n", spacer);
4601 if (SetUpPartitionCounters () == ERROR)
4603 if ((chainParams.stat = (STATS *) SafeCalloc (numTopologies, sizeof (STATS))) == NULL)
4606 memAllocs[ALLOC_STATS] = YES;
4608 /* deal with the reference tree files */
4609 // for (k=0; k<numTopologies; k++)
4610 for (n=0; n<nRefRun; n++)
4613 sprintf (inRefName, "%s.t", comptreeParams.comptFileName2);
4615 sprintf (inRefName, "%s.run%d.t", comptreeParams.comptFileName2, n+1);
4617 /* Examine each ref tree file, save info to tFileInfo */
4618 if (ExamineSumtFile(inRefName, &tFileInfo, treeName, &sumtParams.brlensDef) == ERROR)
4620 if (longestL < tFileInfo.longestLineLength)
4622 longestL = tFileInfo.longestLineLength;
4623 lineBuf = (char *) SafeRealloc (lineBuf, (size_t)longestL * sizeof(char));
4626 MrBayesPrint ("%s Problem allocating string for reading tree file\n", spacer);
4631 /* calculate burnin */
4632 if (chainParams.relativeBurnin == YES)
4633 burnin = (int)(chainParams.burninFraction * tFileInfo.numTreesInLastBlock);
4635 burnin = chainParams.chainBurnIn;
4636 if (burnin >= tFileInfo.numTreesInLastBlock)
4638 MrBayesPrint ("%s Burnin should be smaller than the total number of trees\n", spacer);
4642 /* open the ref tree file */
4643 if ((fpTre = OpenTextFileR (inRefName)) == NULL)
4645 /* ...and fast forward to beginning in last tree block */
4646 for (i=0; i <= tFileInfo.lastTreeBlockBegin; i++)
4648 if (fgets(lineBuf, longestL-2, fpTre) == NULL)
4652 /* process each ref tree */
4653 for (i=1; i <= tFileInfo.numTreesInLastBlock; i++)
4656 if (fgets (lineBuf, longestL-2, fpTre) == NULL)
4658 s = strtok (lineBuf, " ");
4660 while (strcmp (s, "tree") != 0);
4662 /* add reference trees to partition counters, discarding burnin */
4665 s = strtok (NULL, ";");
4670 if (ResetTopology (t, s) == ERROR)
4672 if (AddTreeToPartitionCounters (t, 0, n) == ERROR)
4678 /* close the tree file */
4679 SafeFclose (&fpTre);
4681 /* end reference tree files */
4683 /* open output file */
4684 strcpy (outName, comptreeParams.comptOutfile);
4685 strcat (outName, ".sdsf");
4686 if ((fpOut = OpenNewMBPrintFile (outName)) == NULL)
4688 /* print stamp and header */
4689 if ((int)strlen(stamp) > 1)
4690 MrBayesPrintf (fpOut, "[ID: %s]\n", stamp);
4691 if (chainParams.diagnStat == AVGSTDDEV)
4692 MrBayesPrintf (fpOut, "Gen\tASDSF\n");
4694 MrBayesPrintf (fpOut, "Gen\tMSDSF\n");
4696 /* Examine the tree file to be compared, save info to tFileInfo */
4697 strcpy(inName, comptreeParams.comptFileName1);
4698 if (ExamineSumtFile(inName, &tFileInfo, treeName, &sumtParams.brlensDef) == ERROR)
4700 if (longestL < tFileInfo.longestLineLength)
4702 longestL = tFileInfo.longestLineLength;
4703 lineBuf = (char *) SafeRealloc (lineBuf, (size_t)longestL * sizeof(char));
4706 MrBayesPrint ("%s Problem allocating string for reading tree file\n", spacer);
4711 /* open the tree file to be compared */
4712 if ((fpTre = OpenTextFileR (inName)) == NULL)
4714 /* ...and fast forward to beginning in last tree block */
4715 for (i=0; i <= tFileInfo.lastTreeBlockBegin; i++)
4717 if (fgets(lineBuf, longestL-2, fpTre) == NULL)
4721 /* process each tree to be compared and print SDSF to file */
4722 for (i=1; i <= tFileInfo.numTreesInLastBlock; i++)
4725 if (fgets (lineBuf, longestL-2, fpTre) == NULL)
4727 s = strtok (lineBuf, " ");
4729 while (strcmp (s, "tree") != 0);
4731 s = strtok (NULL, ";");
4732 gen = atoi(s+4); // 4 is offset to get rid of "rep." in tree name
4737 /* add the tree to partition counters */
4738 if (ResetTopology (t, s) == ERROR)
4740 if (AddTreeToPartitionCounters (t, 0, nRefRun) == ERROR)
4744 /* calculate and write stdev of split freq */
4745 CalcTopoConvDiagn2 (nTre);
4746 if (chainParams.stat[0].numPartitions == 0)
4748 MrBayesPrintf (fpOut, "%d\tNA\n", gen);
4750 else if (chainParams.diagnStat == AVGSTDDEV)
4752 MrBayesPrintf (fpOut, "%d\t%lf\n", gen, chainParams.stat[0].avgStdDev);
4756 MrBayesPrintf (fpOut, "%d\t%lf\n", gen, chainParams.stat[0].max);
4760 /* change back to the actual numRuns, end of hack */
4761 chainParams.numRuns -= 1;
4763 /* close tree file */
4764 SafeFclose (&fpTre);
4765 /* close output file */
4766 SafeFclose (&fpOut);
4776 MrBayesPrint ("%s Error in DoCompRefTree\n", spacer);
4778 chainParams.numRuns -= 1;
4779 SafeFclose (&fpTre);
4780 SafeFclose (&fpOut);
4790 #if defined (PRINT_RATEMUL_CPP)
4791 int DELETE_ME_count_taxa(PolyNode *p)
4795 if (p->left==NULL) {
4804 sum+=DELETE_ME_count_taxa(p);
4811 void DELETE_ME_dump_depth(PolyNode *p)
4813 /*print depth of two taxa clade*/
4815 if (p->left != NULL && p->left->left == NULL && p->left->sib != NULL && p->left->sib->left == NULL){
4816 fprintf (rateMultfp,"%f\n",p->depth);
4818 if (p->left != NULL && p->left->left == NULL && p->left->sib != NULL && p->left->sib->left == NULL){
4819 if (p->left->depth > 0.1 && p->left->sib->depth > 0.1)
4820 fprintf (rateMultfp,"%f\n",p->depth);
4823 /* print depth of three taxa clade */
4824 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) ||
4825 (p->left != NULL && p->left->left != NULL && p->left->left->left == NULL && p->left->left->sib->left == NULL && (p->left->sib->left == NULL))) {
4826 if (DELETE_ME_count_taxa(p)==2)
4827 fprintf (rateMultfp,"%f\n",p->depth);
4832 DELETE_ME_dump_depth(p);
4841 int i, j=0, k, n, len, longestName, treeNo, numTreePartsToPrint,
4842 maxWidthID, maxWidthNumberPartitions, maxNumTaxa, tableWidth=0, unreliable, oneUnreliable,
4844 MrBFlt f, var_s, sum_s, stddev_s=0.0, sumsq_s, sumStdDev=0.0, maxStdDev=0.0, sumPSRF=0.0,
4845 maxPSRF=0.0, avgStdDev=0.0, avgPSRF=0.0, min_s=0.0, max_s=0.0, numPSRFSamples=0, min;
4847 char *s=NULL, tempName[120], fileName[120], treeName[100], divString[100];
4848 char *tempStr=NULL; /*not static because error ext is handeled*/
4851 PartCtr **treeParts=NULL,*tmp;
4852 SumtFileInfo sumtFileInfo;
4856 #define SCREEN_WIDTH 80
4858 # if defined (MPI_ENABLED)
4863 /* Ensure that we read trees with sumt code and not user tree code */
4864 inSumtCommand = YES;
4866 /* set file pointers to NULL */
4867 fp = fpParts = fpTstat = fpVstat = fpCon = fpTrees = NULL;
4869 strcpy(treeName,"tree"); //in case if parameter is not specified in a .t file
4871 /* Check if there is anything to do */
4872 if (sumtParams.table == NO && sumtParams.summary == NO && sumtParams.showConsensus == NO)
4874 MrBayesPrint ("%s Nothing to do, all output parameters (Table, Summary, Consensus) set to 'NO'\n", spacer);
4878 SafeStrcpy(&tempStr, " ");
4880 /* Initialize sumtParams struct */
4881 sumtParams.numTaxa = 0;
4882 sumtParams.taxaNames = NULL;
4883 sumtParams.BitsLongsNeeded = 0;
4884 sumtParams.tree = AllocatePolyTree (numTaxa);
4885 sumtParams.nESets = 0;
4886 sumtParams.nBSets = 0;
4887 sumtParams.eSetName = NULL;
4888 sumtParams.bSetName = NULL;
4889 sumtParams.popSizeSet = NO;
4890 sumtParams.popSizeSetName = NULL;
4891 AllocatePolyTreePartitions (sumtParams.tree);
4892 sumtParams.numFileTrees = (int *) SafeCalloc (2*sumtParams.numRuns+2*numTaxa, sizeof(int));
4893 sumtParams.numFileTreesSampled = sumtParams.numFileTrees + sumtParams.numRuns;
4894 sumtParams.order = sumtParams.numFileTrees + 2*sumtParams.numRuns;
4895 sumtParams.absentTaxa = sumtParams.numFileTrees + 2*sumtParams.numRuns + numTaxa;
4896 if (!sumtParams.numFileTrees)
4898 MrBayesPrint ("%s Problems allocating sumtParams.numFileTrees in DoSumt()\n", spacer);
4902 memAllocs[ALLOC_SUMTPARAMS] = YES;
4904 /* Make sure outgroup is set correctly */
4906 for (treeNo = 0; treeNo < sumtParams.numTrees; treeNo++)
4908 /* initialize across-file tree and partition counters */
4909 sumtParams.numTreesSampled = sumtParams.numTreesEncountered = 0;
4910 numUniqueSplitsFound = numUniqueTreesFound = 0;
4912 /* initialize oneUnreliable && unreliable */
4913 oneUnreliable = unreliable = NO;
4915 /* initialize summary statistics */
4922 /* tell user we are ready to go */
4923 if (sumtParams.numTrees > 1)
4924 sprintf (fileName,"%s.tree%d", sumtParams.sumtFileName, treeNo+1);
4926 strcpy (fileName, sumtParams.sumtFileName);
4928 if (sumtParams.numRuns == 1)
4929 MrBayesPrint ("%s Summarizing trees in file \"%s.t\"\n", spacer, fileName);
4930 else if (sumtParams.numRuns == 2)
4931 MrBayesPrint ("%s Summarizing trees in files \"%s.run1.t\" and \"%s.run2.t\"\n", spacer, fileName, fileName);
4932 else if (sumtParams.numRuns > 2)
4933 MrBayesPrint ("%s Summarizing trees in files \"%s.run1.t\", \"%s.run2.t\",...,\"%s.run%d.t\"\n", spacer, fileName, fileName,fileName,sumtParams.numRuns);
4935 if (chainParams.relativeBurnin == YES)
4936 MrBayesPrint ("%s Using relative burnin ('relburnin=yes'), discarding the first %.0f %% of sampled trees\n",
4937 spacer, chainParams.burninFraction*100.0, chainParams.burninFraction);
4939 MrBayesPrint ("%s Using absolute burnin ('relburnin=no'), discarding the first %d sampled trees\n",
4940 spacer, chainParams.chainBurnIn, chainParams.chainBurnIn);
4942 MrBayesPrint ("%s Writing statistics to files %s.<parts|tstat|vstat|trprobs|con>\n", spacer, sumtParams.sumtOutfile);
4944 for (sumtParams.runId=0; sumtParams.runId < sumtParams.numRuns; sumtParams.runId++)
4946 /* initialize tree counters */
4947 sumtParams.numFileTrees[sumtParams.runId] = sumtParams.numFileTreesSampled[sumtParams.runId] = 0;
4949 /* open binary file */
4950 if (sumtParams.numRuns == 1)
4951 sprintf (tempName, "%s.t", fileName);
4953 sprintf (tempName, "%s.run%d.t", fileName, sumtParams.runId+1);
4954 strcpy(sumtParams.curFileName, tempName);
4956 /* tell user we are examining files if for the first run */
4957 if (sumtParams.runId == 0)
4959 if (sumtParams.numRuns > 1 && sumtParams.numTrees > 1)
4960 MrBayesPrint ("%s Examining first file for tree %d ...\n", spacer, treeNo);
4961 else if (sumtParams.numRuns > 1 && sumtParams.numTrees == 1)
4962 MrBayesPrint ("%s Examining first file ...\n", spacer);
4963 else if (sumtParams.numRuns == 1 && sumtParams.numTrees > 1)
4964 MrBayesPrint ("%s Examining file for tree %d ...\n", spacer, treeNo);
4966 MrBayesPrint ("%s Examining file ...\n", spacer);
4970 if (ExamineSumtFile(tempName, &sumtFileInfo, treeName, &sumtParams.brlensDef) == ERROR)
4973 /* capture values */
4974 if (sumtParams.runId == 0)
4976 /* catch lack of sampled trees */
4977 if (chainParams.relativeBurnin == NO && chainParams.chainBurnIn > sumtFileInfo.numTreesInLastBlock)
4979 MrBayesPrint ("%s No trees are sampled as the burnin exceeds the number of trees in last block\n", spacer);
4980 MrBayesPrint ("%s Try setting burnin to a number less than %d\n", spacer, sumtFileInfo.numTreesInLastBlock);
4984 /* tell the user that everything is fine */
4985 if (sumtParams.runId == 0)
4987 if (sumtFileInfo.numTreeBlocks == 1)
4988 MrBayesPrint ("%s Found one tree block in file \"%s\" with %d trees in last block\n",
4989 spacer, tempName, sumtFileInfo.numTreesInLastBlock);
4992 MrBayesPrint ("%s Found %d tree blocks in file \"%s\" with %d trees in last block\n",
4993 spacer, sumtFileInfo.numTreeBlocks, tempName, sumtFileInfo.numTreesInLastBlock);
4994 MrBayesPrint ("%s Only the %d trees in last tree block will be summarized\n", spacer, sumtFileInfo.numTreesInLastBlock);
4996 sumtParams.numTreesInLastBlock = sumtFileInfo.numTreesInLastBlock;
4997 if (sumtParams.numRuns > 1)
4998 MrBayesPrint ("%s Expecting the same number of trees in the last tree block of all files\n", spacer);
4999 if (chainParams.relativeBurnin == NO)
5000 sumtParams.burnin = chainParams.chainBurnIn;
5002 sumtParams.burnin = (int) (sumtFileInfo.numTreesInLastBlock * chainParams.burninFraction);
5006 if (sumtFileInfo.numTreesInLastBlock != sumtParams.numFileTrees[0])
5008 MrBayesPrint ("%s Found %d trees in first file but %d trees in file \"%s\"\n", spacer,
5009 sumtParams.numFileTrees[0],
5010 sumtFileInfo.numTreesInLastBlock,
5016 /* Now we read the file for real. First, allocate a string for reading the file... */
5017 if (sumtParams.runId == 0 && treeNo == 0)
5019 s = (char *) SafeMalloc ((size_t)(sumtFileInfo.longestLineLength) * sizeof(char));
5022 MrBayesPrint ("%s Problem allocating string for reading sumt file\n", spacer);
5029 s = (char *) SafeMalloc ((size_t)(sumtFileInfo.longestLineLength) * sizeof (char));
5032 MrBayesPrint ("%s Problem reallocating string for reading sumt file\n", spacer);
5037 /* ... open the file ... */
5038 if ((fp = OpenTextFileR(tempName)) == NULL)
5041 /* ...and fast forward to beginning of last tree block. */
5042 for (i=0; i<sumtFileInfo.lastTreeBlockBegin+1; i++)
5044 if (fgets (s, sumtFileInfo.longestLineLength-2, fp) == NULL)
5046 printf ("Error in function: %s at line: %d in file: %s", __FUNCTION__, __LINE__, __FILE__);
5050 # if defined (PRINT_RATEMUL_CPP)
5051 sprintf (tempName, "%s.ratemult", chainParams.chainFileName);
5052 if ((rateMultfp=OpenNewMBPrintFile (tempName)) == NULL)
5054 printf ("Error oppening file: %s to write", tempName);
5057 fprintf (rateMultfp,"rateMult_CPP\n");
5060 /* Set up cheap status bar. */
5061 if (sumtParams.runId == 0)
5063 if (sumtParams.numTrees > 1)
5064 MrBayesPrint ("\n%s Tree reading status for tree %d:\n\n", spacer, treeNo+1);
5066 MrBayesPrint ("\n%s Tree reading status:\n\n", spacer);
5067 MrBayesPrint ("%s 0 10 20 30 40 50 60 70 80 90 100\n", spacer);
5068 MrBayesPrint ("%s v-------v-------v-------v-------v-------v-------v-------v-------v-------v-------v\n", spacer);
5069 MrBayesPrint ("%s *", spacer);
5073 /* ...and parse file, tree-by-tree. We are only parsing lines between the "begin trees" and "end" statements.
5074 We don't actually get those lines, however, but rather the lines between those statements. */
5075 expecting = Expecting(COMMAND);
5076 /* We skip the begin trees statement so we need to set up some variables here */
5078 ResetTranslateTable();
5079 for (i=0; i<sumtFileInfo.lastTreeBlockEnd - sumtFileInfo.lastTreeBlockBegin - 1; i++)
5081 if (fgets (s, sumtFileInfo.longestLineLength-2, fp) == NULL)
5083 printf ("Error in function: %s at line: %d in file: %s", __FUNCTION__, __LINE__, __FILE__);
5085 /*MrBayesPrint ("%s", s);*/
5086 if (ParseCommand (s) == ERROR)
5090 ResetTranslateTable();
5092 /* Finish cheap status bar. */
5093 if (sumtParams.runId == sumtParams.numRuns - 1)
5095 if (numAsterices < 80)
5097 for (i=0; i<80 - numAsterices; i++)
5100 MrBayesPrint ("\n\n");
5103 /* print out information on absent and pruned taxa */
5104 if (sumtParams.runId == sumtParams.numRuns - 1 && treeNo == 0)
5105 PrintSumtTaxaInfo ();
5107 /* tell user how many trees were successfully read */
5108 if (sumtParams.numRuns == 1)
5109 MrBayesPrint ("%s Read %d trees from last tree block (sampling %d of them)\n", spacer,
5110 sumtParams.numTreesEncountered, sumtParams.numTreesSampled);
5111 else if (sumtParams.numRuns > 1)
5113 if (sumtParams.runId != 0 && sumtParams.numFileTreesSampled[sumtParams.runId]!=sumtParams.numFileTreesSampled[0])
5115 if (sumtParams.runId == sumtParams.numRuns - 1)
5116 MrBayesPrint ("\n\n");
5117 MrBayesPrint ("%s Found %d post-burnin trees in the first file but %d post-burnin trees in file %d\n",
5119 sumtParams.numFileTreesSampled[0],
5120 sumtParams.numFileTreesSampled[sumtParams.runId],
5121 sumtParams.runId+1);
5124 if (sumtParams.runId == sumtParams.numRuns - 1)
5126 MrBayesPrint ("%s Read a total of %d trees in %d files (sampling %d of them)\n", spacer,
5127 sumtParams.numTreesEncountered,
5128 sumtParams.numRuns, sumtParams.numTreesSampled);
5129 MrBayesPrint ("%s (Each file contained %d trees of which %d were sampled)\n", spacer,
5130 sumtParams.numFileTrees[0],
5131 sumtParams.numFileTreesSampled[0]);
5135 /* Check that at least one tree was read in. */
5136 if (sumtParams.numTreesSampled <= 0)
5138 MrBayesPrint ("%s No trees read in\n", spacer);
5143 } /* next run for this tree */
5145 /* Extract partition counter pointers */
5146 treeParts = (PartCtr **) SafeCalloc ((size_t)numUniqueSplitsFound, sizeof(PartCtr *));
5148 PartCtrUppass(partCtrRoot, treeParts, &i);
5150 min = (sumtParams.minPartFreq * (sumtParams.numTreesSampled/sumtParams.numRuns));
5151 numTreePartsToPrint=numUniqueSplitsFound;
5152 for (i=0; i<numTreePartsToPrint;)
5154 for (j=0; j<sumtParams.numRuns;j++)
5156 if (treeParts[i]->count[j]>=min)
5159 if (j==sumtParams.numRuns)
5161 numTreePartsToPrint--;
5162 tmp=treeParts[numTreePartsToPrint];
5163 treeParts[numTreePartsToPrint]=treeParts[i];
5172 /* Sort taxon partitions (clades, splits) ... */
5173 SortPartCtr (treeParts, 0, numTreePartsToPrint-1);
5175 /* Sort root and tips among those splits always present */
5176 SortTerminalPartCtr (treeParts, numUniqueSplitsFound);
5178 /* open output files for summary information (three files) */
5179 if (OpenSumtFiles (treeNo) == ERROR)
5182 /* Print partitions to screen. */
5186 for (k=0; k<numTaxa; k++)
5188 if (taxaInfo[k].isDeleted == NO && sumtParams.absentTaxa[k] == NO)
5190 len = (int) strlen (taxaNames[k]);
5191 if (len > longestName)
5194 if (sumtParams.table == YES)
5196 MrBayesPrint (" \n");
5197 MrBayesPrint ("%s General explanation: \n", spacer);
5198 MrBayesPrint (" \n");
5199 MrBayesPrint ("%s In an unrooted tree, a taxon bipartition (split) is specified by removing a \n", spacer);
5200 MrBayesPrint ("%s branch, thereby dividing the species into those to the left and those to the \n", spacer);
5201 MrBayesPrint ("%s right of the branch. Here, taxa to one side of the removed branch are denoted \n", spacer);
5202 MrBayesPrint ("%s '.' and those to the other side are denoted '*'. Specifically, the '.' symbol \n", spacer);
5203 MrBayesPrint ("%s is used for the taxa on the same side as the outgroup. \n", spacer);
5204 MrBayesPrint (" \n");
5205 MrBayesPrint ("%s In a rooted or clock tree, the tree is rooted using the model and not by \n", spacer);
5206 MrBayesPrint ("%s reference to an outgroup. Each bipartition therefore corresponds to a clade, \n", spacer);
5207 MrBayesPrint ("%s that is, a group that includes all the descendants of a particular branch in \n", spacer);
5208 MrBayesPrint ("%s the tree. Taxa that are included in each clade are denoted using '*', and \n", spacer);
5209 MrBayesPrint ("%s taxa that are not included are denoted using the '.' symbol. \n", spacer);
5210 MrBayesPrint (" \n");
5211 MrBayesPrint ("%s The output first includes a key to all the bipartitions with frequency larger \n", spacer);
5212 MrBayesPrint ("%s or equual to (Minpartfreq) in at least one run. Minpartfreq is a parameter to \n", spacer);
5213 MrBayesPrint ("%s sumt command and currently it is set to %1.2lf. This is followed by a table \n", spacer, sumtParams.minPartFreq);
5214 MrBayesPrint ("%s with statistics for the informative bipartitions (those including at least \n", spacer);
5215 MrBayesPrint ("%s two taxa), sorted from highest to lowest probability. For each bipartition, \n", spacer);
5216 MrBayesPrint ("%s the table gives the number of times the partition or split was observed in all\n", spacer);
5217 MrBayesPrint ("%s runs (#obs) and the posterior probability of the bipartition (Probab.), which \n", spacer);
5218 MrBayesPrint ("%s is the same as the split frequency. If several runs are summarized, this is \n", spacer);
5219 MrBayesPrint ("%s followed by the minimum split frequency (Min(s)), the maximum frequency \n", spacer);
5220 MrBayesPrint ("%s (Max(s)), and the standard deviation of frequencies (Stddev(s)) across runs. \n", spacer);
5221 MrBayesPrint ("%s The latter value should approach 0 for all bipartitions as MCMC runs converge.\n", spacer);
5222 MrBayesPrint (" \n");
5223 MrBayesPrint ("%s This is followed by a table summarizing branch lengths, node heights (if a \n", spacer);
5224 MrBayesPrint ("%s clock model was used) and relaxed clock parameters (if a relaxed clock model \n", spacer);
5225 MrBayesPrint ("%s was used). The mean, variance, and 95 %% credible interval are given for each \n", spacer);
5226 MrBayesPrint ("%s of these parameters. If several runs are summarized, the potential scale \n", spacer);
5227 MrBayesPrint ("%s reduction factor (PSRF) is also given; it should approach 1 as runs converge. \n", spacer);
5228 MrBayesPrint ("%s Node heights will take calibration points into account, if such points were \n", spacer);
5229 MrBayesPrint ("%s used in the analysis. \n", spacer);
5230 MrBayesPrint ("%s \n", spacer);
5231 MrBayesPrint ("%s Note that Stddev may be unreliable if the partition is not present in all \n", spacer);
5232 MrBayesPrint ("%s runs (the last column indicates the number of runs that sampled the partition \n", spacer);
5233 MrBayesPrint ("%s if more than one run is summarized). The PSRF is not calculated at all if \n", spacer);
5234 MrBayesPrint ("%s the partition is not present in all runs.The PSRF is also sensitive to small \n", spacer);
5235 MrBayesPrint ("%s sample sizes and it should only be considered a rough guide to convergence \n", spacer);
5236 MrBayesPrint ("%s since some of the assumptions allowing one to interpret it as a true potential\n", spacer);
5237 MrBayesPrint ("%s scale reduction factor are violated in MrBayes. \n", spacer);
5238 MrBayesPrint ("%s \n", spacer);
5239 MrBayesPrint ("%s List of taxa in bipartitions: \n", spacer);
5240 MrBayesPrint (" \n");
5242 for (k=0; k<numTaxa; k++)
5244 if (taxaInfo[k].isDeleted == NO && sumtParams.absentTaxa[k] == NO)
5246 MrBayesPrint ("%s %4d -- %s\n", spacer, j++, taxaNames[k]);
5252 if (sumtParams.numTrees > 1 && (sumtParams.table == YES || sumtParams.summary == YES))
5254 MrBayesPrint ("\n\n");
5255 MrBayesPrint ("%s Results for tree number %d\n", spacer, treeNo+1);
5256 MrBayesPrint ("%s ==========================\n\n", spacer);
5259 /* First print key to taxon bipartitions */
5260 if (sumtParams.table == YES)
5262 MrBayesPrint ("\n");
5263 if (sumtParams.numTrees == 1)
5264 MrBayesPrint ("%s Key to taxon bipartitions (saved to file \"%s.parts\"):\n\n", spacer, sumtParams.sumtOutfile);
5266 MrBayesPrint ("%s Key to taxon bipartitions (saved to file \"%s.tree%d.parts\"):\n\n", spacer, sumtParams.sumtOutfile, treeNo+1);
5269 /* calculate a couple of numbers that are handy to have */
5270 /*numTreePartsToPrint = 0;
5271 for (i=0; i<numUniqueSplitsFound; i++)
5273 if ((MrBFlt)treeParts[i]->totCount/(MrBFlt)sumtParams.numTreesSampled >= sumtParams.minPartFreq)
5274 numTreePartsToPrint++;
5277 maxWidthID = (int) (log10 (numTreePartsToPrint)) + 1;
5280 maxNumTaxa = SCREEN_WIDTH - 9;
5282 for (j=0; j<sumtParams.numTaxa; j+=maxNumTaxa)
5284 /* print header to screen and to parts file simultaneously */
5286 /* first print header to screen */
5287 MrBayesPrint ("%s ", spacer);
5289 MrBayesPrint ("%*s -- Partition\n", maxWidthID, "ID");
5291 MrBayesPrint ("%*s -- Partition (continued)\n", maxWidthID, "ID");
5292 tableWidth = maxWidthID + 4 + sumtParams.numTaxa;
5293 if (tableWidth > SCREEN_WIDTH)
5294 tableWidth = SCREEN_WIDTH;
5295 MrBayesPrint ("%s ", spacer);
5296 for (i=0; i<tableWidth; i++)
5298 MrBayesPrint ("\n");
5300 /* now print header to file */
5301 MrBayesPrintf (fpParts, "ID\tPartition\n");
5303 /* now, show partitions that were found on screen; print to .parts file simultaneously */
5304 mask = SafeCalloc (sumtParams.BitsLongsNeeded, sizeof(BitsLong));
5305 for (i=0; i<sumtParams.numTaxa; i++)
5307 for (i=0; i<numTreePartsToPrint; i++)
5310 if (IsBitSet(localOutGroup, x->partition) == YES && sumtParams.isRooted == NO)
5311 FlipBits(x->partition, sumtParams.BitsLongsNeeded, mask);
5313 if ((NumBits(x->partition, sumtParams.BitsLongsNeeded) == numLocalTaxa || NumBits(x->partition, sumtParams.BitsLongsNeeded) == 0) && sumtParams.isClock == NO)
5316 MrBayesPrint ("%s %*d -- ", spacer, maxWidthID, i);
5317 if (sumtParams.numTaxa <= maxNumTaxa)
5318 ShowParts (stdout, x->partition, sumtParams.numTaxa);
5321 if (sumtParams.numTaxa - j > maxNumTaxa)
5322 ShowSomeParts (stdout, x->partition, j, maxNumTaxa);
5324 ShowSomeParts (stdout, x->partition, j, sumtParams.numTaxa - j);
5327 MrBayesPrint ("\n");
5329 MrBayesPrintf (fpParts, "%d\t", i);
5330 ShowParts (fpParts, x->partition, sumtParams.numTaxa);
5331 MrBayesPrintf (fpParts, "\n");
5335 /* finish screen table */
5336 if (sumtParams.table == YES)
5338 MrBayesPrint ("%s ", spacer);
5339 for (i=0; i<tableWidth; i++)
5343 MrBayesPrint ("\n");
5344 if (oneUnreliable == YES)
5346 MrBayesPrint ("%s * The partition was not found in all runs so the values are unreliable\n\n", spacer);
5350 MrBayesPrint ("\n");
5355 /* Second, print statitistics for taxon bipartitions */
5356 if (sumtParams.table == YES)
5358 if (sumtParams.isRooted == NO)
5359 MrBayesPrint ("%s Summary statistics for informative taxon bipartitions\n", spacer);
5361 MrBayesPrint ("%s Summary statistics for informative taxon bipartitions (clades)\n", spacer);
5362 MrBayesPrint ("%s (saved to file \"%s.tstat\"):\n\n", spacer, sumtParams.sumtOutfile);
5365 /* calculate a couple of numbers that are handy to have */
5366 /*numTreePartsToPrint = 0;
5367 for (i=0; i<numUniqueSplitsFound; i++)
5369 if ((MrBFlt)treeParts[i]->totCount/(MrBFlt)sumtParams.numTreesSampled >= sumtParams.minPartFreq)
5370 numTreePartsToPrint++;
5373 maxWidthID = (int) (log10 (numTreePartsToPrint)) + 1;
5376 maxWidthNumberPartitions = (int) (log10 (treeParts[0]->totCount)) + 1;
5377 if (maxWidthNumberPartitions < 4)
5378 maxWidthNumberPartitions = 4;
5380 /* print header to screen and to parts file simultaneously */
5381 if (sumtParams.table == YES)
5383 /* first print header to screen */
5384 MrBayesPrint ("%s ", spacer);
5385 MrBayesPrint ("%*s ", maxWidthID, "ID");
5386 tableWidth = maxWidthID + 3;
5387 MrBayesPrint ("#obs");
5389 for (i=4; i<maxWidthNumberPartitions; i++)
5394 MrBayesPrint (" Probab.");
5396 if (sumtParams.numRuns > 1)
5398 MrBayesPrint (" Sd(s)+ ");
5399 MrBayesPrint (" Min(s) Max(s) ");
5401 MrBayesPrint (" Nruns ");
5404 MrBayesPrint ("\n%s ", spacer);
5405 for (i=0; i<tableWidth; i++)
5409 MrBayesPrint ("\n");
5411 /* now print header to file */
5412 if (sumtParams.numRuns > 1)
5413 MrBayesPrintf (fpTstat, "ID\t#obs\tProbability(=s)\tStddev(s)\tMin(s)\tMax(s)\tNruns\n");
5415 MrBayesPrintf (fpTstat, "ID\t#obs\tProbability(=s)\n");
5418 /* now, show informative partitions that were found on screen; print to .tstat file simultaneously */
5419 for (i=0; i<numTreePartsToPrint; i++)
5423 /* skip uninformative partitions */
5424 if (NumBits(x->partition, sumtParams.BitsLongsNeeded) <= 1 || NumBits(x->partition, sumtParams.BitsLongsNeeded) == sumtParams.numTaxa)
5426 if (NumBits(x->partition, sumtParams.BitsLongsNeeded) == sumtParams.numTaxa - 1 && sumtParams.isRooted == NO)
5429 if (sumtParams.table == YES)
5431 MrBayesPrint ("%s %*d", spacer, maxWidthID, i);
5433 MrBayesPrintf (fpTstat, "%d\t", i);
5435 if (sumtParams.numRuns > 1)
5441 for (n=j=0; n<sumtParams.numRuns; n++)
5443 if (x->count[n] > 0)
5445 f = (MrBFlt) x->count[n] / (MrBFlt) sumtParams.numFileTreesSampled[n];
5453 var_s = sumsq_s - sum_s * sum_s / (MrBFlt) sumtParams.numRuns;
5454 var_s /= (sumtParams.numRuns - 1);
5456 stddev_s = sqrt (var_s);
5459 if (j == sumtParams.numRuns)
5464 oneUnreliable = YES;
5467 if (sumtParams.table == YES)
5469 f = (MrBFlt) x->totCount / (MrBFlt) sumtParams.numTreesSampled;
5470 MrBayesPrint (" %*d %1.6lf", maxWidthNumberPartitions, x->totCount, f);
5471 MrBayesPrintf (fpTstat, "\t%d\t%s", x->totCount, MbPrintNum(f));
5472 if (sumtParams.numRuns > 1)
5474 MrBayesPrint (" %1.6lf %1.6lf %1.6lf", stddev_s, min_s, max_s);
5475 MrBayesPrint (" %3d", j);
5476 MrBayesPrintf (fpTstat, "\t%s", MbPrintNum(stddev_s));
5477 MrBayesPrintf (fpTstat, "\t%s", MbPrintNum(min_s));
5478 MrBayesPrintf (fpTstat, "\t%s", MbPrintNum(max_s));
5479 MrBayesPrintf (fpTstat, "\t%d", j);
5481 MrBayesPrintf (fpTstat, "\n");
5482 if (unreliable == YES)
5483 MrBayesPrint (" *\n");
5485 MrBayesPrint ("\n");
5486 sumStdDev += stddev_s;
5487 if (stddev_s > maxStdDev)
5488 maxStdDev = stddev_s;
5492 /* finish screen table */
5493 if (sumtParams.table == YES)
5495 MrBayesPrint ("%s ", spacer);
5496 for (i=0; i<tableWidth; i++)
5500 MrBayesPrint ("\n");
5501 if (sumtParams.numRuns > 1)
5503 MrBayesPrint ("%s + Convergence diagnostic (standard deviation of split frequencies)\n", spacer);
5504 MrBayesPrint ("%s should approach 0.0 as runs converge.\n\n", spacer);
5506 if (oneUnreliable == YES)
5507 MrBayesPrint ("%s * The partition was not found in all runs so the values are unreliable\n", spacer);
5510 /* Third, print statitistics for branch and node parameters */
5511 if (sumtParams.table == YES)
5513 MrBayesPrint ("\n");
5514 MrBayesPrint ("%s Summary statistics for branch and node parameters\n", spacer);
5515 MrBayesPrint ("%s (saved to file \"%s.vstat\"):\n", spacer, sumtParams.sumtOutfile);
5518 if (sumtParams.table == YES)
5520 /* calculate longest header */
5521 longestHeader = 9; /* length of 'parameter' */
5522 i = (int)(log10(numTreePartsToPrint)) + 3; /* length of partition specifier including [] */
5523 len = i + (int)(strlen(treeName)) + 2; /* length of length{m}[n] or height{m}[n] */
5524 if (len > longestHeader)
5525 longestHeader = len;
5526 for (j=0; j<sumtParams.nBSets; j++)
5528 len = (int) strlen(sumtParams.tree->bSetName[j]) + 7 + i;
5529 if (len > longestHeader)
5530 longestHeader = len;
5532 for (j=0; j<sumtParams.nESets; j++)
5534 len = (int) strlen(sumtParams.tree->eSetName[j]) + 8 + i;
5535 if (len > longestHeader)
5536 longestHeader = len;
5538 if (sumtParams.popSizeSet == YES)
5540 len = (int) strlen(sumtParams.tree->popSizeSetName) + i;
5541 if (len > longestHeader)
5542 longestHeader = len;
5545 /* print the header rows */
5546 MrBayesPrint ("\n");
5547 if (sumtParams.HPD == NO)
5548 MrBayesPrint ("%s %*c 95%% Cred. Interval\n", spacer, longestHeader, ' ');
5550 MrBayesPrint ("%s %*c 95%% HPD Interval\n", spacer, longestHeader, ' ');
5551 MrBayesPrint ("%s %*c --------------------\n", spacer, longestHeader, ' ');
5553 MrBayesPrint ("%s Parameter%*c Mean Variance Lower Upper Median", spacer, longestHeader-9, ' ');
5554 tableWidth = 68 + longestHeader - 9;
5555 if (sumtParams.HPD == YES)
5556 MrBayesPrintf (fpVstat, "Parameter\tMean\tVariance\tCredInt_Lower\tCredInt_Upper\tMedian", spacer, longestHeader-9, ' ');
5558 MrBayesPrintf (fpVstat, "Parameter\tMean\tVariance\tHPD_Lower\tHPD_Upper\tMedian", spacer, longestHeader-9, ' ');
5559 if (sumtParams.numRuns > 1)
5561 MrBayesPrint (" PSRF+ Nruns");
5563 MrBayesPrintf (fpVstat, "\tPSRF\tNruns");
5565 MrBayesPrint ("\n");
5566 MrBayesPrintf (fpVstat, "\n");
5568 MrBayesPrint ("%s ", spacer);
5569 for (j=0; j<tableWidth; j++)
5573 MrBayesPrint ("\n");
5576 strcpy (divString, treeName+4);
5577 for (i=1; i<numTreePartsToPrint; i++)
5580 tempStrLength=(int)strlen(tempStr);
5581 SafeSprintf (&tempStr,&tempStrLength, "length%s[%d]", divString, i);
5582 len = (int) strlen(tempStr);
5584 GetSummary (x->length, sumtParams.numRuns, x->count, &theStats, sumtParams.HPD);
5586 MrBayesPrint ("%s %-*s ", spacer, longestHeader, tempStr);
5587 MrBayesPrintf (fpVstat, "%s", tempStr);
5589 PrintSumtTableLine(sumtParams.numRuns, x->count, &theStats, &numPSRFSamples, &maxPSRF, &sumPSRF);
5594 if (sumtParams.isClock == YES)
5596 strcpy (divString, treeName+4);
5597 for (i=0; i<numTreePartsToPrint; i++)
5600 tempStrLength=(int)strlen(tempStr);
5601 SafeSprintf (&tempStr,&tempStrLength, "height%s[%d]", divString, i);
5602 len = (int) strlen(tempStr);
5604 GetSummary (x->height, sumtParams.numRuns, x->count, &theStats, sumtParams.HPD);
5606 MrBayesPrint ("%s %-*s ", spacer, longestHeader, tempStr);
5607 MrBayesPrintf (fpVstat, "%s", tempStr);
5609 PrintSumtTableLine(sumtParams.numRuns, x->count, &theStats, &numPSRFSamples, &maxPSRF, &sumPSRF);
5614 if (sumtParams.isCalibrated == YES)
5616 strcpy (divString, treeName+4);
5617 for (i=0; i<numTreePartsToPrint; i++)
5620 tempStrLength=(int)strlen(tempStr);
5621 SafeSprintf (&tempStr,&tempStrLength, "age%s[%d]", divString, i);
5622 len = (int) strlen(tempStr);
5624 GetSummary (x->age, sumtParams.numRuns, x->count, &theStats, sumtParams.HPD);
5626 MrBayesPrint ("%s %-*s ", spacer, longestHeader, tempStr);
5627 MrBayesPrintf (fpVstat, "%s", tempStr);
5629 PrintSumtTableLine(sumtParams.numRuns, x->count, &theStats, &numPSRFSamples, &maxPSRF, &sumPSRF);
5633 /* print effective branch lengths */
5634 if (sumtParams.isRelaxed == YES)
5636 for (i=0; i<sumtParams.nBSets; i++)
5638 for (j=1; j<numTreePartsToPrint; j++)
5641 tempStrLength=(int)strlen(tempStr);
5642 SafeSprintf (&tempStr,&tempStrLength, "%s_length[%d]", sumtParams.bSetName[i], j);
5643 len = (int) strlen(tempStr);
5645 GetSummary (x->bLen[i], sumtParams.numRuns, x->count, &theStats, sumtParams.HPD);
5647 MrBayesPrint ("%s %-*s ", spacer, longestHeader, tempStr);
5648 MrBayesPrintf (fpVstat, "%s", tempStr);
5650 PrintSumtTableLine(sumtParams.numRuns, x->count, &theStats, &numPSRFSamples, &maxPSRF, &sumPSRF);
5652 for (j=1; j<numTreePartsToPrint; j++)
5655 tempStrLength=(int)strlen(tempStr);
5656 SafeSprintf (&tempStr,&tempStrLength, "%s_rate[%d]", sumtParams.bSetName[i], j);
5657 len = (int) strlen(tempStr);
5659 GetSummary (x->bRate[i], sumtParams.numRuns, x->count, &theStats, sumtParams.HPD);
5661 MrBayesPrint ("%s %-*s ", spacer, longestHeader, tempStr);
5662 MrBayesPrintf (fpVstat, "%s", tempStr);
5664 PrintSumtTableLine(sumtParams.numRuns, x->count, &theStats, &numPSRFSamples, &maxPSRF, &sumPSRF);
5667 for (i=0; i<sumtParams.nESets; i++)
5669 for (j=1; j<numTreePartsToPrint; j++)
5672 tempStrLength=(int)strlen(tempStr);
5673 SafeSprintf (&tempStr,&tempStrLength, "%s_nEvents[%d]", sumtParams.eSetName[i], j);
5674 len = (int) strlen(tempStr);
5676 GetIntSummary (x->nEvents[i], sumtParams.numRuns, x->count, &theStats, sumtParams.HPD);
5678 MrBayesPrint ("%s %-*s ", spacer, longestHeader, tempStr);
5679 MrBayesPrintf (fpVstat, "%s", tempStr);
5681 PrintSumtTableLine(sumtParams.numRuns, x->count, &theStats, &numPSRFSamples, &maxPSRF, &sumPSRF);
5686 /* print population size sets */
5687 if (sumtParams.popSizeSet == YES)
5689 for (j=1; j<numTreePartsToPrint; j++)
5692 tempStrLength=(int)strlen(tempStr);
5693 SafeSprintf (&tempStr,&tempStrLength, "%s[%d]", sumtParams.popSizeSetName, j);
5694 len = (int) strlen(tempStr);
5696 GetSummary (x->popSize, sumtParams.numRuns, x->count, &theStats, sumtParams.HPD);
5698 MrBayesPrint ("%s %-*s ", spacer, longestHeader, tempStr);
5699 MrBayesPrintf (fpVstat, "%s", tempStr);
5701 PrintSumtTableLine(sumtParams.numRuns, x->count, &theStats, &numPSRFSamples, &maxPSRF, &sumPSRF);
5706 MrBayesPrint ("%s ", spacer);
5707 for (j=0; j<tableWidth; j++)
5709 MrBayesPrint ("\n");
5711 if (sumtParams.numRuns > 1)
5713 MrBayesPrint ("%s + Convergence diagnostic (PSRF = Potential Scale Reduction Factor; Gelman\n", spacer);
5714 MrBayesPrint ("%s and Rubin, 1992) should approach 1.0 as runs converge. NA is reported when\n", spacer);
5715 MrBayesPrint ("%s deviation of parameter values within all runs is 0 or when a parameter\n", spacer);
5716 MrBayesPrint ("%s value (a branch length, for instance) is not sampled in all runs.\n", spacer);
5719 if (oneUnreliable == YES)
5721 MrBayesPrint ("%s * The partition was not found in all runs so the values are unreliable.\n", spacer);
5723 MrBayesPrint ("\n\n");
5726 /* Exclude trivial splits when calculating average standard deviation of split frequencies. */
5727 avgStdDev = sumStdDev / (numTreePartsToPrint-sumtParams.numTaxa-1);
5728 avgPSRF = sumPSRF / numPSRFSamples;
5730 if (sumtParams.numRuns > 1 && sumtParams.summary == YES)
5732 MrBayesPrint ("%s Summary statistics for partitions with frequency >= %1.2lf in at least one run:\n", spacer, sumtParams.minPartFreq);
5733 MrBayesPrint ("%s Average standard deviation of split frequencies = %1.6lf\n", spacer, avgStdDev);
5734 MrBayesPrint ("%s Maximum standard deviation of split frequencies = %1.6lf\n", spacer, maxStdDev);
5736 if (sumtParams.brlensDef == YES && sumtParams.numRuns > 1 && sumtParams.summary == YES)
5738 MrBayesPrint ("%s Average PSRF for parameter values (excluding NA and >10.0) = %1.3lf\n", spacer, avgPSRF);
5740 MrBayesPrint ("%s Maximum PSRF for parameter values = NA\n", spacer);
5742 MrBayesPrint ("%s Maximum PSRF for parameter values = %1.3lf\n", spacer, maxPSRF);
5744 MrBayesPrint ("\n");
5746 SortPartCtr (treeParts, 0, numUniqueSplitsFound-1); /* We sort again but this time we sort all partitions instead of just first numTreePartsToPrintNow */
5747 /* make the majority rule consensus tree */
5748 if (sumtParams.showConsensus == YES && ConTree (treeParts, numUniqueSplitsFound) == ERROR)
5751 /* get probabilities of individual trees */
5752 if (TreeProb () == ERROR)
5756 if (sumtParams.printBrlensToFile == YES && PrintBrlensToFile (treeParts, numUniqueSplitsFound, treeNo) == ERROR)
5760 SafeFclose (&fpParts);
5761 SafeFclose (&fpTstat);
5762 SafeFclose (&fpVstat);
5763 SafeFclose (&fpCon);
5764 SafeFclose (&fpTrees);
5766 # if defined (PRINT_RATEMUL_CPP)
5767 SafeFclose (&rateMultfp);
5770 /* free pointer array to partitions */
5773 FreePartCtr (partCtrRoot);
5775 FreeTreeCtr (treeCtrRoot);
5779 /* free memory and file pointers */
5783 /* reset numLocalTaxa and localOutGroup */
5786 # if defined (MPI_ENABLED)
5790 expecting = Expecting(COMMAND);
5792 SafeFree ((void **)&tempStr);
5798 /* free sumtParams */
5802 /* close files in case they are open*/
5804 SafeFclose (&fpParts);
5805 SafeFclose (&fpTstat);
5806 SafeFclose (&fpVstat);
5807 SafeFclose (&fpCon);
5808 SafeFclose (&fpTrees);
5810 # if defined (PRINT_RATEMUL_CPP)
5811 SafeFclose (&rateMultfp);
5814 /* free pointer array to partitions, part and tree counters */
5816 FreePartCtr (partCtrRoot);
5817 FreeTreeCtr (treeCtrRoot);
5821 /* reset taxon set */
5824 expecting = Expecting(COMMAND);
5826 SafeFree ((void **)&tempStr);
5832 int DoSumtParm (char *parmName, char *tkn)
5838 if (defMatrix == NO)
5840 MrBayesPrint ("%s A matrix must be specified before sumt can be used\n", spacer);
5844 if (expecting == Expecting(PARAMETER))
5846 expecting = Expecting(EQUALSIGN);
5850 if (!strcmp(parmName, "Xxxxxxxxxx"))
5852 expecting = Expecting(PARAMETER);
5853 expecting |= Expecting(SEMICOLON);
5855 /* set Filename (sumtParams.sumtFileName) ***************************************************/
5856 else if (!strcmp(parmName, "Filename"))
5858 if (expecting == Expecting(EQUALSIGN))
5860 expecting = Expecting(ALPHA);
5863 else if (expecting == Expecting(ALPHA))
5867 MrBayesPrint ("%s Maximum allowed length of file name is 99 characters. The given name:\n", spacer);
5868 MrBayesPrint ("%s '%s'\n", spacer,tkn);
5869 MrBayesPrint ("%s has %d characters.\n", spacer,strlen(tkn));
5872 strcpy (sumtParams.sumtFileName, tkn);
5873 strcpy(sumtParams.sumtOutfile, tkn);
5874 MrBayesPrint ("%s Setting sumt filename and outputname to %s\n", spacer, sumtParams.sumtFileName);
5875 expecting = Expecting(PARAMETER) | Expecting(SEMICOLON);
5880 /* set Relburnin (chainParams.relativeBurnin) ********************************************************/
5881 else if (!strcmp(parmName, "Relburnin"))
5883 if (expecting == Expecting(EQUALSIGN))
5884 expecting = Expecting(ALPHA);
5885 else if (expecting == Expecting(ALPHA))
5887 if (IsArgValid(tkn, tempStr) == NO_ERROR)
5889 if (!strcmp(tempStr, "Yes"))
5890 chainParams.relativeBurnin = YES;
5892 chainParams.relativeBurnin = NO;
5896 MrBayesPrint ("%s Invalid argument for Relburnin\n", spacer);
5899 if (chainParams.relativeBurnin == YES)
5900 MrBayesPrint ("%s Using relative burnin (a fraction of samples discarded).\n", spacer);
5902 MrBayesPrint ("%s Using absolute burnin (a fixed number of samples discarded).\n", spacer);
5903 expecting = Expecting(PARAMETER) | Expecting(SEMICOLON);
5910 /* set Burnin (chainParams.chainBurnIn) ***********************************************************/
5911 else if (!strcmp(parmName, "Burnin"))
5913 if (expecting == Expecting(EQUALSIGN))
5914 expecting = Expecting(NUMBER);
5915 else if (expecting == Expecting(NUMBER))
5917 sscanf (tkn, "%d", &tempI);
5918 chainParams.chainBurnIn = tempI;
5919 MrBayesPrint ("%s Setting urn-in to %d\n", spacer, chainParams.chainBurnIn);
5920 expecting = Expecting(PARAMETER) | Expecting(SEMICOLON);
5927 /* set Burninfrac (chainParams.burninFraction) ************************************************************/
5928 else if (!strcmp(parmName, "Burninfrac"))
5930 if (expecting == Expecting(EQUALSIGN))
5931 expecting = Expecting(NUMBER);
5932 else if (expecting == Expecting(NUMBER))
5934 sscanf (tkn, "%lf", &tempD);
5937 MrBayesPrint ("%s Burnin fraction too low (< 0.01)\n", spacer);
5942 MrBayesPrint ("%s Burnin fraction too high (> 0.50)\n", spacer);
5945 chainParams.burninFraction = tempD;
5946 MrBayesPrint ("%s Setting burnin fraction to %.2f\n", spacer, chainParams.burninFraction);
5947 expecting = Expecting(PARAMETER) | Expecting(SEMICOLON);
5954 /* set Nruns (sumtParams.numRuns) *******************************************************/
5955 else if (!strcmp(parmName, "Nruns"))
5957 if (expecting == Expecting(EQUALSIGN))
5958 expecting = Expecting(NUMBER);
5959 else if (expecting == Expecting(NUMBER))
5961 sscanf (tkn, "%d", &tempI);
5964 MrBayesPrint ("%s Nruns must be at least 1\n", spacer);
5969 sumtParams.numRuns = tempI;
5970 MrBayesPrint ("%s Setting sumt nruns to %d\n", spacer, sumtParams.numRuns);
5971 expecting = Expecting(PARAMETER) | Expecting(SEMICOLON);
5977 /* set Ntrees (sumtParams.numTrees) *******************************************************/
5978 else if (!strcmp(parmName, "Ntrees"))
5980 if (expecting == Expecting(EQUALSIGN))
5981 expecting = Expecting(NUMBER);
5982 else if (expecting == Expecting(NUMBER))
5984 sscanf (tkn, "%d", &tempI);
5987 MrBayesPrint ("%s Ntrees must be at least 1\n", spacer);
5992 sumtParams.numTrees = tempI;
5993 MrBayesPrint ("%s Setting sumt ntrees to %d\n", spacer, sumtParams.numTrees);
5994 expecting = Expecting(PARAMETER) | Expecting(SEMICOLON);
6000 /* set Contype (sumtParams.sumtConType) *****************************************************/
6001 else if (!strcmp(parmName, "Contype"))
6003 if (expecting == Expecting(EQUALSIGN))
6004 expecting = Expecting(ALPHA);
6005 else if (expecting == Expecting(ALPHA))
6007 if (IsArgValid(tkn, tempStr) == NO_ERROR)
6009 strcpy (sumtParams.sumtConType, tempStr);
6010 MrBayesPrint ("%s Setting sumt contype to %s\n", spacer, sumtParams.sumtConType);
6012 expecting = Expecting(PARAMETER) | Expecting(SEMICOLON);
6017 /* set Conformat (sumtParams.consensusFormat) *****************************************************/
6018 else if (!strcmp(parmName, "Conformat"))
6020 if (expecting == Expecting(EQUALSIGN))
6021 expecting = Expecting(ALPHA);
6022 else if (expecting == Expecting(ALPHA))
6024 if (IsArgValid(tkn, tempStr) == NO_ERROR)
6026 if (!strcmp(tempStr,"Figtree"))
6028 MrBayesPrint ("%s Setting sumt conformat to Figtree\n", spacer);
6029 sumtParams.consensusFormat = FIGTREE;
6033 MrBayesPrint ("%s Setting sumt conformat to Simple\n", spacer);
6034 sumtParams.consensusFormat = SIMPLE;
6039 MrBayesPrint ("%s Invalid argument for calctreeprobs\n", spacer);
6042 expecting = Expecting(PARAMETER) | Expecting(SEMICOLON);
6047 /* set Calctreeprobs (sumtParams.calcTreeprobs) *********************************************/
6048 else if (!strcmp(parmName, "Calctreeprobs"))
6050 if (expecting == Expecting(EQUALSIGN))
6051 expecting = Expecting(ALPHA);
6052 else if (expecting == Expecting(ALPHA))
6054 if (IsArgValid(tkn, tempStr) == NO_ERROR)
6056 if (!strcmp(tempStr, "Yes"))
6057 sumtParams.calcTreeprobs = YES;
6059 sumtParams.calcTreeprobs = NO;
6063 MrBayesPrint ("%s Invalid argument for calctreeprobs\n", spacer);
6066 if (sumtParams.calcTreeprobs == YES)
6067 MrBayesPrint ("%s Setting calctreeprobs to yes\n", spacer);
6069 MrBayesPrint ("%s Setting calctreeprobs to no\n", spacer);
6070 expecting = Expecting(PARAMETER) | Expecting(SEMICOLON);
6075 /* set Showtreeprobs (sumtParams.showSumtTrees) *********************************************/
6076 else if (!strcmp(parmName, "Showtreeprobs"))
6078 if (expecting == Expecting(EQUALSIGN))
6079 expecting = Expecting(ALPHA);
6080 else if (expecting == Expecting(ALPHA))
6082 if (IsArgValid(tkn, tempStr) == NO_ERROR)
6084 if (!strcmp(tempStr, "Yes"))
6085 sumtParams.showSumtTrees = YES;
6087 sumtParams.showSumtTrees = NO;
6091 MrBayesPrint ("%s Invalid argument for showtreeprobs\n", spacer);
6094 if (sumtParams.showSumtTrees == YES)
6095 MrBayesPrint ("%s Setting showtreeprobs to yes\n", spacer);
6097 MrBayesPrint ("%s Setting showtreeprobs to no\n", spacer);
6098 expecting = Expecting(PARAMETER) | Expecting(SEMICOLON);
6103 /* set Hpd (sumpParams.HPD) ********************************************************/
6104 else if (!strcmp(parmName, "Hpd"))
6106 if (expecting == Expecting(EQUALSIGN))
6107 expecting = Expecting(ALPHA);
6108 else if (expecting == Expecting(ALPHA))
6110 if (IsArgValid(tkn, tempStr) == NO_ERROR)
6112 if (!strcmp(tempStr, "Yes"))
6113 sumtParams.HPD = YES;
6115 sumtParams.HPD = NO;
6119 MrBayesPrint ("%s Invalid argument for Hpd\n", spacer);
6122 if (sumtParams.HPD == YES)
6123 MrBayesPrint ("%s Reporting 95 %% region of Highest Posterior Density (HPD).\n", spacer);
6125 MrBayesPrint ("%s Reporting median interval containing 95 %% of values.\n", spacer);
6126 expecting = Expecting(PARAMETER) | Expecting(SEMICOLON);
6131 /* set Printbrlens (sumtParams.printBrlensToFile) *********************************************/
6132 else if (!strcmp(parmName, "Printbrlens"))
6134 if (expecting == Expecting(EQUALSIGN))
6135 expecting = Expecting(ALPHA);
6136 else if (expecting == Expecting(ALPHA))
6138 if (IsArgValid(tkn, tempStr) == NO_ERROR)
6140 if (!strcmp(tempStr, "Yes"))
6141 sumtParams.printBrlensToFile = YES;
6143 sumtParams.printBrlensToFile = NO;
6147 MrBayesPrint ("%s Invalid argument for printbrlens\n", spacer);
6150 if (sumtParams.printBrlensToFile == YES)
6151 MrBayesPrint ("%s Setting printbrlens to yes\n", spacer);
6153 MrBayesPrint ("%s Setting printbrlens to no\n", spacer);
6154 expecting = Expecting(PARAMETER) | Expecting(SEMICOLON);
6159 /* set Brlensgeq (sumtParams.brlensFreqDisplay) *******************************************************/
6160 else if (!strcmp(parmName, "Brlensgeq"))
6162 if (expecting == Expecting(EQUALSIGN))
6163 expecting = Expecting(NUMBER);
6164 else if (expecting == Expecting(NUMBER))
6166 sscanf (tkn, "%lf", &tempD);
6167 sumtParams.brlensFreqDisplay = tempD;
6168 MrBayesPrint ("%s Printing branch lengths to file for partitions with probability >= %lf\n", spacer, sumtParams.brlensFreqDisplay);
6169 expecting = Expecting(PARAMETER) | Expecting(SEMICOLON);
6174 /* set Ordertaxa (sumtParams.orderTaxa) *********************************************/
6175 else if (!strcmp(parmName, "Ordertaxa"))
6177 if (expecting == Expecting(EQUALSIGN))
6178 expecting = Expecting(ALPHA);
6179 else if (expecting == Expecting(ALPHA))
6181 if (IsArgValid(tkn, tempStr) == NO_ERROR)
6183 if (!strcmp(tempStr, "Yes"))
6184 sumtParams.orderTaxa = YES;
6186 sumtParams.orderTaxa = NO;
6190 MrBayesPrint ("%s Invalid argument for ordertaxa\n", spacer);
6193 if (sumtParams.orderTaxa == YES)
6194 MrBayesPrint ("%s Setting ordertaxa to yes\n", spacer);
6196 MrBayesPrint ("%s Setting ordertaxa to no\n", spacer);
6197 expecting = Expecting(PARAMETER) | Expecting(SEMICOLON);
6202 /* set Outputname (sumtParams.sumtOutfile) *******************************************************/
6203 else if (!strcmp(parmName, "Outputname"))
6205 if (expecting == Expecting(EQUALSIGN))
6207 expecting = Expecting(ALPHA);
6210 else if (expecting == Expecting(ALPHA))
6212 sscanf (tkn, "%s", tempStr);
6213 strcpy (sumtParams.sumtOutfile, tempStr);
6214 MrBayesPrint ("%s Setting sumt output file name to \"%s\"\n", spacer, sumtParams.sumtOutfile);
6215 expecting = Expecting(PARAMETER) | Expecting(SEMICOLON);
6220 /* set Table (sumtParams.table) ********************************************************/
6221 else if (!strcmp(parmName, "Table"))
6223 if (expecting == Expecting(EQUALSIGN))
6224 expecting = Expecting(ALPHA);
6225 else if (expecting == Expecting(ALPHA))
6227 if (IsArgValid(tkn, tempStr) == NO_ERROR)
6229 if (!strcmp(tempStr, "Yes"))
6230 sumtParams.table = YES;
6232 sumtParams.table = NO;
6236 MrBayesPrint ("%s Invalid argument for Table (valid arguments are 'yes' and 'no')\n", spacer);
6239 if (sumtParams.table == YES)
6240 MrBayesPrint ("%s Setting sumt to compute table of partition frequencies\n", spacer);
6242 MrBayesPrint ("%s Setting sumt not to compute table of partition frequencies\n", spacer);
6243 expecting = Expecting(PARAMETER) | Expecting(SEMICOLON);
6248 /* set Summary (sumtParams.summary) ********************************************************/
6249 else if (!strcmp(parmName, "Summary"))
6251 if (expecting == Expecting(EQUALSIGN))
6252 expecting = Expecting(ALPHA);
6253 else if (expecting == Expecting(ALPHA))
6255 if (IsArgValid(tkn, tempStr) == NO_ERROR)
6257 if (!strcmp(tempStr, "Yes"))
6258 sumtParams.summary = YES;
6260 sumtParams.summary = NO;
6264 MrBayesPrint ("%s Invalid argument for 'Summary' (valid arguments are 'yes' and 'no')\n", spacer);
6267 if (sumtParams.summary == YES)
6268 MrBayesPrint ("%s Setting sumt to summary statistics\n", spacer);
6270 MrBayesPrint ("%s Setting sumt not to compute summary statistics\n", spacer);
6271 expecting = Expecting(PARAMETER) | Expecting(SEMICOLON);
6276 /* set Consensus (sumtParams.showConsensus) ********************************************************/
6277 else if (!strcmp(parmName, "Consensus"))
6279 if (expecting == Expecting(EQUALSIGN))
6280 expecting = Expecting(ALPHA);
6281 else if (expecting == Expecting(ALPHA))
6283 if (IsArgValid(tkn, tempStr) == NO_ERROR)
6285 if (!strcmp(tempStr, "Yes"))
6286 sumtParams.showConsensus = YES;
6288 sumtParams.showConsensus = NO;
6292 MrBayesPrint ("%s Invalid argument for Consensus (valid arguments are 'yes' and 'no')\n", spacer);
6295 if (sumtParams.showConsensus == YES)
6296 MrBayesPrint ("%s Setting sumt to show consensus trees\n", spacer);
6298 MrBayesPrint ("%s Setting sumt not to show consensus trees\n", spacer);
6299 expecting = Expecting(PARAMETER) | Expecting(SEMICOLON);
6304 /* set Minpartfreq (sumtParams.minPartFreq) *******************************************************/
6305 else if (!strcmp(parmName, "Minpartfreq"))
6307 if (expecting == Expecting(EQUALSIGN))
6308 expecting = Expecting(NUMBER);
6309 else if (expecting == Expecting(NUMBER))
6311 sscanf (tkn, "%lf", &tempD);
6312 sumtParams.minPartFreq = tempD;
6313 MrBayesPrint ("%s Including partitions with probability greater than or equal to %lf in summary statistics\n", spacer, sumtParams.minPartFreq);
6314 expecting = Expecting(PARAMETER) | Expecting(SEMICOLON);
6327 int DoSumtTree (void)
6329 int i, j, z, printEvery, nAstPerPrint, burnin;
6334 # if defined (PRINT_RATEMUL_CPP)
6335 /* get depths if relevant */
6336 if (sumtParams.tree->isClock)
6337 GetPolyDepths (sumtParams.tree);
6338 if (rateMultfp !=NULL && sumtParams.tree->root !=NULL)
6339 DELETE_ME_dump_depth(sumtParams.tree->root);
6342 /* increment number of trees read in */
6343 sumtParams.numFileTrees[sumtParams.runId]++;
6344 sumtParams.numTreesEncountered++;
6346 /* update status bar */
6347 if (sumtParams.numTreesInLastBlock * sumtParams.numRuns < 80)
6350 nAstPerPrint = 80 / (sumtParams.numTreesInLastBlock * sumtParams.numRuns);
6351 if (sumtParams.numTreesEncountered % printEvery == 0)
6353 for (i=0; i<nAstPerPrint; i++)
6362 x = (MrBFlt)(sumtParams.numTreesInLastBlock * sumtParams.numRuns) / (MrBFlt) (80);
6363 y = (MrBFlt)(sumtParams.numFileTrees[sumtParams.runId] + sumtParams.numTreesInLastBlock * sumtParams.runId) / x;
6365 if (numAsterices < z)
6373 if (inComparetreeCommand == YES)
6374 burnin = comptreeParams.burnin;
6376 burnin = sumtParams.burnin;
6378 if (sumtParams.numFileTrees[sumtParams.runId] > burnin)
6380 /* increment the number of trees sampled */
6381 sumtParams.numFileTreesSampled[sumtParams.runId]++;
6382 sumtParams.numTreesSampled++;
6384 /* get the tree we just read in */
6385 t = sumtParams.tree;
6387 /* move calculation root for nonrooted trees if necessary */
6388 MovePolyCalculationRoot (t, localOutGroup);
6390 /* check taxon set and outgroup */
6391 if (sumtParams.runId == 0 && sumtParams.numFileTreesSampled[0] == 1)
6393 if (isTranslateDef == YES && isTranslateDiff == YES)
6395 /* we are using a translate block with different taxa set */
6396 if (t->nNodes - t->nIntNodes != numTranslates)
6398 MrBayesPrint ("%s ERROR: Expected %d taxa; found %d taxa\n", spacer, numTranslates, sumtParams.numTaxa);
6401 for (i=0; i<numTaxa; i++)
6402 sumtParams.absentTaxa[i] = NO;
6403 for (i=numTranslates; i<numTaxa; i++)
6404 sumtParams.absentTaxa[i] = YES;
6405 for (i=0; i<numTranslates; i++)
6406 if (strcmp(transFrom[i], taxaNames[outGroupNum]) == 0)
6408 if (i == numTranslates)
6409 localOutGroup = 0; /* no previous outgroup assignment is valid */
6415 /* we are using the current taxa set */
6416 for (i=0; i<numTaxa; i++)
6417 sumtParams.absentTaxa[i] = YES;
6418 for (i=0; i<t->nNodes; i++)
6420 p = t->allDownPass[i];
6421 if (p->left == NULL)
6422 sumtParams.absentTaxa[p->index] = NO;
6425 for (i=j=0; i<numTaxa; i++)
6427 if (sumtParams.absentTaxa[i] == NO && taxaInfo[i].isDeleted == NO)
6429 if (i == outGroupNum)
6436 /* now we can safely prune the tree based on taxaInfo[].isDeleted */
6437 /* note that PrunePolyTree relies on labels to recognize tips */
6440 /* reset tip and int node indices in case some taxa deleted */
6441 ResetTipIndices (t);
6442 ResetIntNodeIndices(t);
6444 /* set basic parameters */
6445 sumtParams.numTaxa = t->nNodes - t->nIntNodes;
6446 numLocalTaxa = sumtParams.numTaxa;
6447 sumtParams.BitsLongsNeeded = ((numLocalTaxa-1) / nBitsInALong) + 1;
6448 if (t->isRooted == YES)
6449 sumtParams.orderLen = numLocalTaxa - 2;
6451 sumtParams.orderLen = numLocalTaxa - 3;
6455 /* the following block was conditioned with if (isTranslateDef == NO || isTranslateDiff == NO)
6456 The reason was not clearly stated but it prevents exclusion of taxa to work in case when the condition does not hold.
6457 My guess is that before PrunePolyTree() relied on indeses of tips be set as in original matrix.
6458 Now it is not needed after PrunePolyTree and ResetTipIndices ware modified to use labels istead of indexes to recognize tips.*/
6460 for (i=0; i<t->nNodes; i++)
6462 p = t->allDownPass[i];
6463 if (p->left == NULL && taxaInfo[p->index].isDeleted == NO && sumtParams.absentTaxa[p->index] == YES)
6465 MrBayesPrint ("%s Taxon %d should not be in sampled tree\n", spacer, p->index + 1);
6470 /* now we can safely prune the tree based on taxaInfo[].isDeleted */
6473 /* reset tip and int node indices in case some taxa deleted */
6474 ResetTipIndices (t);
6475 ResetIntNodeIndices(t);
6478 /* check that all taxa are included */
6479 if (t->nNodes - t->nIntNodes != sumtParams.numTaxa)
6481 MrBayesPrint ("%s Expecting %d taxa but tree '%s' in file '%s' has %d taxa\n",
6482 spacer, sumtParams.numTaxa, t->name, sumtParams.curFileName, t->nNodes-t->nIntNodes);
6487 if (sumtParams.runId == 0 && sumtParams.numFileTreesSampled[0] == 1)
6489 /* harvest labels (can only be done safely after pruning) */
6490 for (i=0; i<sumtParams.numTaxa; i++)
6492 for (j=0; j<t->nNodes; j++)
6494 p = t->allDownPass[j];
6495 if (p->index == i) {
6496 if (strlen(p->label)>99)
6498 MrBayesPrint ("%s Taxon name %s is too long. Maximun 99 characters is allowed.\n", spacer, p->label);
6501 AddString(&sumtParams.taxaNames, i, p->label);
6507 /* check that tree agrees with template */
6508 if (sumtParams.numTreesSampled == 1)
6510 sumtParams.brlensDef = t->brlensDef;
6511 sumtParams.isRooted = t->isRooted;
6512 sumtParams.isClock = t->isClock;
6513 sumtParams.isCalibrated = t->isCalibrated;
6514 sumtParams.isRelaxed = t->isRelaxed;
6515 sumtParams.nBSets = 0;
6516 sumtParams.nESets = 0;
6517 for (i=0; i<t->nBSets; i++)
6518 AddString(&sumtParams.bSetName,sumtParams.nBSets++,t->bSetName[i]);
6519 for (i=0; i<t->nESets; i++)
6520 AddString(&sumtParams.eSetName,sumtParams.nESets++,t->eSetName[i]);
6521 if (t->popSizeSet == YES)
6523 sumtParams.popSizeSet = YES;
6524 sumtParams.popSizeSetName = (char *) SafeCalloc (strlen(t->popSizeSetName)+1, sizeof(char));
6525 strcpy(sumtParams.popSizeSetName, t->popSizeSetName);
6528 sumtParams.popSizeSet = NO;
6530 else /* if (sumtParams.numTreesSampled > 1) */
6532 if (sumtParams.brlensDef != t->brlensDef)
6534 MrBayesPrint ("%s Trees with and without branch lengths mixed\n", spacer);
6537 if (sumtParams.isRooted != t->isRooted)
6539 if (sumtParams.isRooted == YES)
6540 MrBayesPrint ("%s Expected rooted tree but tree '%s' in file '%s' is not rooted\n",
6541 spacer, t->name, sumtParams.curFileName);
6542 else if (sumtParams.isRooted == NO)
6543 MrBayesPrint ("%s Expected unrooted tree but tree '%s' in file '%s' is rooted\n",
6544 spacer, t->name, sumtParams.curFileName);
6547 if (sumtParams.isClock != t->isClock)
6549 if (sumtParams.isClock == YES)
6550 MrBayesPrint ("%s Expected clock tree but tree '%s' in file '%s' is not clock\n",
6551 spacer, t->name, sumtParams.curFileName);
6552 else if (sumtParams.isClock == NO)
6553 MrBayesPrint ("%s Expected nonclock tree but tree '%s' in file '%s' is clock\n",
6554 spacer, t->name, sumtParams.curFileName);
6557 if (sumtParams.isCalibrated != t->isCalibrated)
6559 if (sumtParams.isCalibrated == YES)
6560 MrBayesPrint ("%s Expected calibrated tree but tree '%s' in file '%s' is not calibrated\n",
6561 spacer, t->name, sumtParams.curFileName);
6562 else if (sumtParams.isCalibrated == NO)
6563 MrBayesPrint ("%s Expected noncalibrated tree but tree '%s' in file '%s' is calibrated\n",
6564 spacer, t->name, sumtParams.curFileName);
6567 if (inComparetreeCommand == NO && sumtParams.isRelaxed != t->isRelaxed)
6569 if (sumtParams.isRelaxed == YES)
6570 MrBayesPrint ("%s Expected relaxed clock tree but tree '%s' in file '%s' is not relaxed\n",
6571 spacer, t->name, sumtParams.curFileName);
6572 else if (sumtParams.isRelaxed == NO)
6573 MrBayesPrint ("%s Expected unrooted tree but tree '%s' in file '%s' is rooted\n",
6574 spacer, t->name, sumtParams.curFileName);
6577 if (inComparetreeCommand == NO && (sumtParams.nESets != t->nESets || sumtParams.nBSets != t->nBSets))
6579 MrBayesPrint ("%s Tree '%s' in file '%s' does not have the expected relaxed clock parameters\n",
6580 spacer, t->name, sumtParams.curFileName);
6585 /* set partitions for tree */
6586 ResetPolyTreePartitions(t);
6588 /* get depths if relevant */
6592 /* get ages if relevant */
6593 if (t->isCalibrated)
6596 /* add partitions to counters */
6597 for (i=0; i<t->nNodes; i++)
6599 p = t->allDownPass[i];
6600 partCtrRoot = AddSumtPartition (partCtrRoot, t, p, sumtParams.runId);
6603 /* add the tree to relevant tree list */
6604 if (inSumtCommand == YES)
6606 if (t->isRooted == YES)
6607 StoreRPolyTopology (t, sumtParams.order);
6608 else /* if (sumtParams.isRooted == NO) */
6609 StoreUPolyTopology (t, sumtParams.order);
6610 treeCtrRoot = AddSumtTree (treeCtrRoot, sumtParams.order);
6614 i = sumtParams.numFileTreesSampled[sumtParams.runId] - 1;
6615 if (StoreSumtTree (packedTreeList[sumtParams.runId], i, t) == ERROR)
6619 /* Display the tree nodes. */
6629 int ExamineSumtFile (char *fileName, SumtFileInfo *sumtFileInfo, char *treeName, int *brlensDef)
6631 int i, foundBegin, lineTerm, inTreeBlock, blockErrors, inSumtComment, lineNum, numTreesInBlock,
6633 char sumtToken[100], *s, *sumtTokenP;
6636 /* open binary file */
6637 if ((fp = OpenBinaryFileR(fileName)) == NULL)
6640 /* find out what type of line termination is used for file 1 */
6641 lineTerm = LineTermType (fp);
6642 if (lineTerm != LINETERM_MAC && lineTerm != LINETERM_DOS && lineTerm != LINETERM_UNIX)
6644 MrBayesPrint ("%s Unknown line termination for file \"%s\"\n", spacer, fileName);
6648 /* find length of longest line in either file */
6649 sumtFileInfo->longestLineLength = LongestLine (fp);
6650 sumtFileInfo->longestLineLength += 10;
6652 /* allocate a string long enough to hold a line */
6653 s = (char *)SafeMalloc((size_t)(sumtFileInfo->longestLineLength) * sizeof(char));
6656 MrBayesPrint ("%s Problem allocating string for examining file \"%s\"\n", spacer, fileName);
6660 /* close binary file */
6663 foundBegin = inTreeBlock = blockErrors = inSumtComment = NO;
6664 lineNum = numTreesInBlock = 0;
6665 sumtFileInfo->numTreeBlocks = 0;
6666 sumtFileInfo->lastTreeBlockBegin = 0;
6667 sumtFileInfo->lastTreeBlockEnd = 0;
6668 sumtFileInfo->numTreesInLastBlock = 0;
6670 /* open text file */
6671 if ((fp = OpenTextFileR(fileName))==NULL)
6673 MrBayesPrint ("%s Could not read file \"%s\" in text mode \n", spacer, fileName);
6678 while (fgets (s, sumtFileInfo->longestLineLength-2, fp) != NULL)
6683 if (GetToken (sumtToken, &tokenType, &sumtTokenP))
6685 if (IsSame("[", sumtToken) == SAME)
6686 inSumtComment = YES;
6687 if (IsSame("]", sumtToken) == SAME)
6690 if (inSumtComment == YES)
6692 if (IsSame ("Param", sumtToken) == SAME)
6694 /* extract the tree name */
6695 if (GetToken (sumtToken, &tokenType, &sumtTokenP)) /* get the colon */
6697 if (GetToken (sumtToken, &tokenType, &sumtTokenP)) /* get the tree name */
6699 strcpy (treeName, sumtToken);
6700 if (GetToken (sumtToken, &tokenType, &sumtTokenP))
6702 while (IsSame("]", sumtToken) != SAME)
6704 strcat (treeName, sumtToken);
6705 if (GetToken (sumtToken, &tokenType, &sumtTokenP))
6711 else /* if (inSumtComment == NO) */
6713 if (foundBegin == YES)
6715 if (IsSame("Trees", sumtToken) == SAME)
6717 numTreesInBlock = 0;
6720 sumtFileInfo->lastTreeBlockBegin = lineNum;
6725 if (IsSame("Begin", sumtToken) == SAME)
6727 if (foundBegin == YES)
6729 MrBayesPrint ("%s Found inappropriate \"Begin\" statement in file\n", spacer);
6734 else if (IsSame("End", sumtToken) == SAME)
6736 if (inTreeBlock == YES)
6738 sumtFileInfo->numTreeBlocks++;
6740 sumtFileInfo->lastTreeBlockEnd = lineNum;
6744 MrBayesPrint ("%s Found inappropriate \"End\" statement in file\n", spacer);
6747 sumtFileInfo->numTreesInLastBlock = numTreesInBlock;
6749 else if (IsSame("Tree", sumtToken) == SAME)
6751 if (inTreeBlock == YES)
6754 if (numTreesInBlock == 1)
6757 for (i=0; s[i]!='\0'; i++)
6769 MrBayesPrint ("%s Found a \"Tree\" statement that is not in a tree block\n", spacer);
6776 } while (*sumtToken);
6780 /* Now, check some aspects of the tree file, such as the number of tree blocks and whether they are properly terminated. */
6781 if (inTreeBlock == YES)
6783 MrBayesPrint ("%s Unterminated tree block in file %s. You probably need to\n", spacer, fileName);
6784 MrBayesPrint ("%s add a new line to the end of the file with \"End;\" on it.\n", spacer);
6787 if (inSumtComment == YES)
6789 MrBayesPrint ("%s Unterminated comment in file %s\n", spacer, fileName);
6792 if (blockErrors == YES)
6794 MrBayesPrint ("%s Found formatting errors in file %s\n", spacer, fileName);
6797 if (sumtFileInfo->lastTreeBlockEnd < sumtFileInfo->lastTreeBlockBegin)
6799 MrBayesPrint ("%s Problem reading tree file %s\n", spacer, fileName);
6802 if (sumtFileInfo->numTreesInLastBlock <= 0)
6804 MrBayesPrint ("%s No trees were found in last tree block of file %s\n", spacer, fileName);
6816 /* FreePartCtr: Recursively free partition counter nodes */
6817 void FreePartCtr (PartCtr *r)
6824 FreePartCtr (r->left);
6825 FreePartCtr (r->right);
6827 /* free relaxed clock parameters: eRate, nEvents, bRate */
6828 if (sumtParams.nESets > 0)
6830 for (i=0; i<sumtParams.nESets; i++)
6832 for (j=0; j<sumtParams.numRuns; j++)
6833 free (r->nEvents[i][j]);
6834 free (r->nEvents[i]);
6838 if (sumtParams.nBSets > 0)
6840 for (i=0; i<sumtParams.nBSets; i++)
6842 for (j=0; j<sumtParams.numRuns; j++)
6844 free (r->bLen [i][j]);
6845 free (r->bRate[i][j]);
6854 /* free basic parameters */
6855 for (i=0; i<sumtParams.numRuns; i++)
6856 free (r->length[i]);
6860 free (r->partition);
6862 numUniqueSplitsFound--;
6867 /* FreeSumtParams: Free parameters allocated in sumtParams struct */
6868 void FreeSumtParams(void)
6872 if (memAllocs[ALLOC_SUMTPARAMS] == YES)
6874 for (i=0; i<sumtParams.numTaxa; i++)
6875 free(sumtParams.taxaNames[i]);
6876 free (sumtParams.taxaNames);
6877 sumtParams.taxaNames = NULL;
6878 if (sumtParams.numFileTrees) free (sumtParams.numFileTrees);
6879 sumtParams.numFileTrees = NULL;
6880 FreePolyTree (sumtParams.tree);
6881 sumtParams.tree = NULL;
6882 if (sumtParams.nBSets > 0)
6884 for (i=0; i<sumtParams.nBSets; i++)
6885 free(sumtParams.bSetName[i]);
6886 free (sumtParams.bSetName);
6887 sumtParams.bSetName = NULL;
6888 sumtParams.nBSets = 0;
6890 if (sumtParams.nESets > 0)
6892 for (i=0; i<sumtParams.nESets; i++)
6893 free(sumtParams.eSetName[i]);
6894 free (sumtParams.eSetName);
6895 sumtParams.eSetName = NULL;
6896 sumtParams.nESets = 0;
6898 if (sumtParams.popSizeSet == YES)
6900 free (sumtParams.popSizeSetName);
6901 sumtParams.popSizeSetName = NULL;
6902 sumtParams.popSizeSet = NO;
6904 memAllocs[ALLOC_SUMTPARAMS] = NO;
6909 /* FreeTreeCtr: Recursively free tree counter nodes */
6910 void FreeTreeCtr (TreeCtr *r)
6915 FreeTreeCtr (r->left);
6916 FreeTreeCtr (r->right);
6920 numUniqueTreesFound--;
6925 /* Label: Calculate length of label and fill in char *label if not NULL */
6926 int Label (PolyNode *p, int addIndex, char *label, int maxLength)
6928 int i, j0, j1, k, n, length, nameLength, index;
6933 /* first calculate length */
6934 if (inSumtCommand == YES && isTranslateDiff == NO)
6936 for (index=i=0; index<numTaxa; index++)
6938 if (sumtParams.absentTaxa[index] == YES || taxaInfo[index].isDeleted == YES)
6950 length = (int)(strlen(p->label)) + 4 + (int)(log10(index+1));
6952 length = (int)(strlen(p->label));
6953 length = (length > maxLength ? maxLength : length);
6955 /* fill in label if label != NULL */
6959 nameLength = length - 4 - (int)(log10(index+1));
6961 nameLength = length;
6963 for (i=0; i<nameLength-1; i++)
6964 label[i] = p->label[i];
6965 if ((int)strlen(p->label) > nameLength)
6968 label[i] = p->label[i];
6975 k = (int)(log10(n)) + 1;
6978 j0 = (int)(log10(n));
6979 j1 = (int)(pow(10,j0));
6980 label[++i] = '0' + n/j1;
6998 int OpenComptFiles (void)
7000 int len, previousFiles, oldNoWarn, oldAutoOverwrite;
7001 char pFilename[120], dFilename[120];
7005 oldAutoOverwrite = autoOverwrite;
7007 /* set file names */
7008 strcpy (pFilename, comptreeParams.comptOutfile);
7009 strcpy (dFilename, comptreeParams.comptOutfile);
7010 strcat (pFilename, ".pairs");
7011 strcat (dFilename, ".dists");
7013 /* one overwrite check for both files */
7017 if ((fpTemp = OpenTextFileR(pFilename)) != NULL)
7019 previousFiles = YES;
7022 if ((fpTemp = OpenTextFileR(dFilename)) != NULL)
7024 previousFiles = YES;
7027 if (previousFiles == YES)
7029 MrBayesPrint("%s There are previous compare results saved using the same filenames.\n", spacer);
7030 if (WantTo("Do you want to overwrite these results") == YES)
7034 autoOverwrite = YES;
7039 MrBayesPrint("%s Please specify a different output file name before running the comparetree command.\n", spacer);
7040 MrBayesPrint("%s You can do that using 'comparetree outputfile=<name>'. You can also move or\n", spacer);
7041 MrBayesPrint("%s rename the old result files.\n", spacer);
7047 if ((fpParts = OpenNewMBPrintFile (pFilename)) == NULL)
7050 autoOverwrite = oldAutoOverwrite;
7053 if ((fpDists = OpenNewMBPrintFile (dFilename)) == NULL)
7056 autoOverwrite = oldAutoOverwrite;
7060 /* Reset file flags */
7062 autoOverwrite = oldAutoOverwrite;
7064 /* print unique identifiers to each file */
7065 len = (int) strlen (stamp);
7068 MrBayesPrintf (fpParts, "[ID: %s]\n", stamp);
7069 MrBayesPrintf (fpDists, "[ID: %s]\n", stamp);
7076 int OpenSumtFiles (int treeNo)
7078 int i, len, oldNoWarn, oldAutoOverwrite, previousFiles;
7079 char pFilename[120], sFilename[120], vFilename[120], cFilename[120], tFilename[120];
7083 oldAutoOverwrite = autoOverwrite;
7085 /* one overwrite check for all files */
7086 if (noWarn == NO && treeNo == 0)
7089 for (i=0; i<sumtParams.numTrees; i++)
7091 if (sumtParams.numTrees > 1)
7093 sprintf (pFilename, "%s.tree%d.parts", sumtParams.sumtOutfile, i+1);
7094 sprintf (sFilename, "%s.tree%d.tstat", sumtParams.sumtOutfile, i+1);
7095 sprintf (vFilename, "%s.tree%d.vstat", sumtParams.sumtOutfile, i+1);
7096 sprintf (cFilename, "%s.tree%d.con.tre", sumtParams.sumtOutfile, i+1);
7097 sprintf (tFilename, "%s.tree%d.trprobs", sumtParams.sumtOutfile, i+1);
7101 sprintf (pFilename, "%s.parts", sumtParams.sumtOutfile);
7102 sprintf (sFilename, "%s.tstat", sumtParams.sumtOutfile);
7103 sprintf (vFilename, "%s.vstat", sumtParams.sumtOutfile);
7104 sprintf (cFilename, "%s.con.tre", sumtParams.sumtOutfile);
7105 sprintf (tFilename, "%s.trprobs", sumtParams.sumtOutfile);
7107 if ((fpTemp = TestOpenTextFileR(pFilename)) != NULL)
7109 previousFiles = YES;
7112 if ((fpTemp = TestOpenTextFileR(sFilename)) != NULL)
7114 previousFiles = YES;
7117 if ((fpTemp = TestOpenTextFileR(vFilename)) != NULL)
7119 previousFiles = YES;
7122 if ((fpTemp = TestOpenTextFileR(cFilename)) != NULL)
7124 previousFiles = YES;
7127 if ((fpTemp = TestOpenTextFileR(tFilename)) != NULL)
7129 previousFiles = YES;
7132 if (previousFiles == YES)
7135 MrBayesPrint("%s There are previous tree sample summaries saved using the same filenames.\n", spacer);
7136 if (WantTo("Do you want to overwrite these results") == YES)
7140 autoOverwrite = YES;
7145 MrBayesPrint("%s Please specify a different output file name before running the sumt command.\n", spacer);
7146 MrBayesPrint("%s You can do that using 'sumt outputfile=<name>'. You can also move or\n", spacer);
7147 MrBayesPrint("%s rename the old result files.\n", spacer);
7154 /* set file names */
7155 if (sumtParams.numTrees > 1)
7157 sprintf (pFilename, "%s.tree%d.parts", sumtParams.sumtOutfile, treeNo+1);
7158 sprintf (sFilename, "%s.tree%d.tstat", sumtParams.sumtOutfile, treeNo+1);
7159 sprintf (vFilename, "%s.tree%d.vstat", sumtParams.sumtOutfile, treeNo+1);
7160 sprintf (cFilename, "%s.tree%d.con.tre", sumtParams.sumtOutfile, treeNo+1);
7161 sprintf (tFilename, "%s.tree%d.trprobs", sumtParams.sumtOutfile, treeNo+1);
7165 sprintf (pFilename, "%s.parts", sumtParams.sumtOutfile);
7166 sprintf (sFilename, "%s.tstat", sumtParams.sumtOutfile);
7167 sprintf (vFilename, "%s.vstat", sumtParams.sumtOutfile);
7168 sprintf (cFilename, "%s.con.tre", sumtParams.sumtOutfile);
7169 sprintf (tFilename, "%s.trprobs", sumtParams.sumtOutfile);
7172 /* open files checking for over-write as appropriate */
7173 if ((fpParts = OpenNewMBPrintFile(pFilename)) == NULL)
7175 if ((fpTstat = OpenNewMBPrintFile(sFilename)) == NULL)
7177 SafeFclose (&fpParts);
7180 if ((fpVstat = OpenNewMBPrintFile(vFilename)) == NULL)
7182 SafeFclose (&fpParts);
7183 SafeFclose (&fpTstat);
7186 if ((fpCon = OpenNewMBPrintFile(cFilename)) == NULL)
7188 SafeFclose (&fpParts);
7189 SafeFclose (&fpTstat);
7190 SafeFclose (&fpVstat);
7193 if (sumtParams.calcTreeprobs == YES)
7195 if ((fpTrees = OpenNewMBPrintFile(tFilename)) == NULL)
7197 SafeFclose (&fpParts);
7198 SafeFclose (&fpTstat);
7199 SafeFclose (&fpVstat);
7200 SafeFclose (&fpCon);
7205 /* print #NEXUS if appropriate */
7206 MrBayesPrintf (fpCon, "#NEXUS\n");
7207 if (sumtParams.calcTreeprobs == YES)
7208 MrBayesPrintf (fpTrees, "#NEXUS\n");
7210 /* print unique identifiers to each file */
7211 len = (int) strlen (stamp);
7214 MrBayesPrintf (fpParts, "[ID: %s]\n", stamp);
7215 MrBayesPrintf (fpTstat, "[ID: %s]\n", stamp);
7216 MrBayesPrintf (fpVstat, "[ID: %s]\n", stamp);
7217 MrBayesPrintf (fpCon, "[ID: %s]\n", stamp);
7218 if (sumtParams.calcTreeprobs == YES)
7219 MrBayesPrintf (fpTrees, "[ID: %s]\n", stamp);
7222 /* Reset noWarn and autoOverwrite */
7223 if (treeNo == sumtParams.numTrees - 1)
7226 autoOverwrite = oldAutoOverwrite;
7233 void PartCtrUppass (PartCtr *r, PartCtr **uppass, int *index)
7237 uppass[(*index)++] = r;
7239 PartCtrUppass (r->left, uppass, index);
7240 PartCtrUppass (r->right, uppass, index);
7245 /* PrintBrlensToFile: Print brlens to file */
7246 int PrintBrlensToFile (PartCtr **treeParts, int numTreeParts, int treeNo)
7248 int i, j, runNo, numBrlens;
7254 if (sumtParams.numTrees > 1)
7255 sprintf (filename, "%s.tree%d.brlens", sumtParams.sumtOutfile, treeNo+1);
7257 sprintf (filename, "%s.brlens", sumtParams.sumtOutfile);
7259 /* Open file checking for over-write as appropriate */
7260 if ((fp = OpenNewMBPrintFile(filename)) == NULL)
7263 /* count number of brlens to print */
7264 for (i=0; i<numTreeParts; i++)
7266 if (treeParts[i]->totCount < sumtParams.brlensFreqDisplay)
7272 for (i=0; i<numBrlens; i++)
7274 MrBayesPrintf (fp, "v[%d]", i+1);
7276 MrBayesPrintf (fp, "\n");
7278 MrBayesPrintf (fp, "\t");
7282 for (i=0; i<numBrlens; i++)
7284 x = treeParts[numBrlens];
7285 for (runNo=0; runNo<sumtParams.numRuns; runNo++)
7287 MrBayesPrintf (fp, "%s", MbPrintNum (x->length[runNo][0]));
7288 for (j=1; j<x->count[i]; j++)
7290 MrBayesPrintf (fp, "\t%s", MbPrintNum (x->length[runNo][j]));
7293 MrBayesPrintf (fp, "\n");
7300 /* PrintConTree: Print consensus tree in standard format readable by TreeView, Paup etc */
7301 void PrintConTree (FILE *fp, PolyTree *t)
7303 MrBayesPrintf (fp, " [Note: This tree contains information on the topology, \n");
7304 MrBayesPrintf (fp, " branch lengths (if present), and the probability\n");
7305 MrBayesPrintf (fp, " of the partition indicated by the branch.]\n");
7306 if (!strcmp(sumtParams.sumtConType, "Halfcompat"))
7307 MrBayesPrintf (fp, " tree con_50_majrule = ");
7309 MrBayesPrintf (fp, " tree con_all_compat = ");
7310 WriteConTree (t->root, fp, YES);
7311 MrBayesPrintf (fp, ";\n");
7312 if (sumtParams.brlensDef == YES)
7314 MrBayesPrintf (fp, "\n");
7315 MrBayesPrintf (fp, " [Note: This tree contains information only on the topology\n");
7316 MrBayesPrintf (fp, " and branch lengths (median of the posterior probability density).]\n");
7317 if (!strcmp(sumtParams.sumtConType, "Halfcompat"))
7318 MrBayesPrintf (fp, " tree con_50_majrule = ");
7320 MrBayesPrintf (fp, " tree con_all_compat = ");
7321 WriteConTree (t->root, fp, NO);
7322 MrBayesPrintf (fp, ";\n");
7327 /* PrintFigTreeConTree: Print consensus tree in rich format for FigTree */
7328 void PrintFigTreeConTree (FILE *fp, PolyTree *t, PartCtr **treeParts)
7330 if (!strcmp(sumtParams.sumtConType, "Halfcompat"))
7331 MrBayesPrintf (fp, " tree con_50_majrule = ");
7333 MrBayesPrintf (fp, " tree con_all_compat = ");
7334 if (t->isRooted == YES)
7335 MrBayesPrintf (fp, "[&R] ");
7337 MrBayesPrintf (fp, "[&U] ");
7339 WriteFigTreeConTree (t->root, fp, treeParts);
7340 MrBayesPrintf (fp, ";\n");
7344 void PrintFigTreeNodeInfo (FILE *fp, PartCtr *x, MrBFlt length)
7346 int i, postProbPercent, postProbSdPercent;
7347 MrBFlt *support, mean, var, min, max;
7350 support = SafeCalloc (sumtParams.numRuns, sizeof(MrBFlt));
7351 for (i=0; i<sumtParams.numRuns; i++)
7353 support[i] = (MrBFlt) x->count[i] / (MrBFlt) sumtParams.numFileTreesSampled[i];
7355 if (sumtParams.numRuns > 1)
7357 MeanVariance (support, sumtParams.numRuns, &mean, &var);
7358 Range (support, sumtParams.numRuns, &min, &max);
7359 postProbPercent = (int) (100.0*mean + 0.5);
7360 postProbSdPercent = (int) (100.0 * sqrt(var) + 0.5);
7361 fprintf (fp, "[&prob=%.8le,prob_stddev=%.8le,prob_range={%.8le,%.8le},prob(percent)=\"%d\",prob+-sd=\"%d+-%d\"",
7362 mean, sqrt(var), min, max, postProbPercent, postProbPercent, postProbSdPercent);
7366 postProbPercent = (int) (100.0*support[0] + 0.5);
7367 fprintf (fp, "[&prob=%.8le,prob(percent)=\"%d\"", support[0], postProbPercent);
7369 if (sumtParams.isClock == YES)
7371 GetSummary (x->height, sumtParams.numRuns, x->count, &theStats, sumtParams.HPD);
7372 if (sumtParams.HPD == YES)
7373 fprintf (fp, ",height_mean=%.8le,height_median=%.8le,height_95%%HPD={%.8le,%.8le}", theStats.mean, theStats.median, theStats.lower, theStats.upper);
7375 fprintf (fp, ",height_mean=%.8le,height_median=%.8le,height_95%%CredInt={%.8le,%.8le}", theStats.mean, theStats.median, theStats.lower, theStats.upper);
7377 if (sumtParams.isCalibrated == YES)
7379 GetSummary (x->age, sumtParams.numRuns, x->count, &theStats, sumtParams.HPD);
7380 if (sumtParams.HPD == YES)
7381 fprintf (fp, ",age_mean=%.8le,age_median=%.8le,age_95%%HPD={%.8le,%.8le}", theStats.mean, theStats.median, theStats.lower, theStats.upper);
7383 fprintf (fp, ",age_mean=%.8le,age_median=%.8le,age_95%%CredInt={%.8le,%.8le}", theStats.mean, theStats.median, theStats.lower, theStats.upper);
7387 fprintf (fp, ":%s", MbPrintNum(length));
7388 if (sumtParams.brlensDef == YES)
7390 GetSummary (x->length, sumtParams.numRuns, x->count, &theStats, sumtParams.HPD);
7391 if (sumtParams.HPD == YES)
7392 fprintf (fp, "[&length_mean=%.8le,length_median=%.8le,length_95%%HPD={%.8le,%.8le}", theStats.mean, theStats.median, theStats.lower, theStats.upper);
7394 fprintf (fp, "[&length_mean=%.8le,length_median=%.8le,length_95%%CredInt={%.8le,%.8le}", theStats.mean, theStats.median, theStats.lower, theStats.upper);
7396 if (sumtParams.isClock == YES && sumtParams.isRelaxed == YES)
7398 for (i=0; i<sumtParams.nBSets; i++)
7400 GetSummary (x->bLen[i], sumtParams.numRuns, x->count, &theStats, sumtParams.HPD);
7401 if (sumtParams.HPD == YES)
7402 fprintf (fp, ",effectivebrlen%s_mean=%lf,effectivebrlen%s_median=%lf,effectivebrlen%s_95%%HPD={%lf,%lf}",
7403 sumtParams.tree->bSetName[i], theStats.mean,
7404 sumtParams.tree->bSetName[i], theStats.median,
7405 sumtParams.tree->bSetName[i], theStats.lower,
7408 fprintf (fp, ",effectivebrlen%s_mean=%lf,effectivebrlen%s_median=%lf,effectivebrlen%s_95%%CredInt={%lf,%lf}",
7409 sumtParams.tree->bSetName[i], theStats.mean,
7410 sumtParams.tree->bSetName[i], theStats.median,
7411 sumtParams.tree->bSetName[i], theStats.lower,
7413 GetSummary (x->bRate[i], sumtParams.numRuns, x->count, &theStats, sumtParams.HPD);
7414 if (sumtParams.HPD == YES)
7415 fprintf (fp, ",rate%s_mean=%lf,rate%s_median=%lf,rate%s_95%%HPD={%lf,%lf}",
7416 sumtParams.tree->bSetName[i], theStats.mean,
7417 sumtParams.tree->bSetName[i], theStats.median,
7418 sumtParams.tree->bSetName[i], theStats.lower,
7421 fprintf (fp, ",rate%s_mean=%lf,rate%s_median=%lf,rate%s_95%%CredInt={%lf,%lf}",
7422 sumtParams.tree->bSetName[i], theStats.mean,
7423 sumtParams.tree->bSetName[i], theStats.median,
7424 sumtParams.tree->bSetName[i], theStats.lower,
7427 for (i=0; i<sumtParams.nESets; i++)
7429 GetIntSummary (x->nEvents[i], sumtParams.numRuns, x->count, &theStats, sumtParams.HPD);
7430 if (sumtParams.HPD == YES)
7431 fprintf (fp, ",nEvents%s_mean=%lf,nEvents%s_median=%lf,nEvents%s_95%%HPD={%lf,%lf}",
7432 sumtParams.tree->eSetName[i], theStats.mean,
7433 sumtParams.tree->eSetName[i], theStats.median,
7434 sumtParams.tree->eSetName[i], theStats.lower,
7437 fprintf (fp, ",nEvents%s_mean=%lf,nEvents%s_median=%lf,nEvents%s_95%%CredInt={%lf,%lf}",
7438 sumtParams.tree->eSetName[i], theStats.mean,
7439 sumtParams.tree->eSetName[i], theStats.median,
7440 sumtParams.tree->eSetName[i], theStats.lower,
7444 if (sumtParams.brlensDef == YES)
7451 void PrintSumtTableLine(int numRuns, int *rowCount, Stat *theStats, MrBFlt *numPSRFSamples, MrBFlt *maxPSRF, MrBFlt *sumPSRF)
7455 MrBayesPrint ("%10.6lf %10.6lf %10.6lf %10.6lf %10.6lf", theStats->mean, theStats->var, theStats->lower, theStats->upper, theStats->median);
7457 MrBayesPrintf (fpVstat, "\t%s", MbPrintNum(theStats->mean));
7458 MrBayesPrintf (fpVstat, "\t%s", MbPrintNum(theStats->var));
7459 MrBayesPrintf (fpVstat, "\t%s", MbPrintNum(theStats->lower));
7460 MrBayesPrintf (fpVstat, "\t%s", MbPrintNum(theStats->upper));
7461 MrBayesPrintf (fpVstat, "\t%s", MbPrintNum(theStats->median));
7465 for (j=k=0; j<numRuns; j++)
7466 if (rowCount[j] > 0)
7468 if (theStats->PSRF < 0.0)
7470 MrBayesPrint (" NA %3d", k);
7471 MrBayesPrintf (fpVstat, "\tNA\t%d", k);
7475 if (theStats->PSRF > 10.0)
7477 MrBayesPrint (" >10.0 %3d", k);
7478 MrBayesPrintf (fpVstat, "\tNA\t%d", k);
7483 MrBayesPrint (" %7.3lf %3d", theStats->PSRF, k);
7484 MrBayesPrintf (fpVstat, "\t%s\t%d", MbPrintNum(theStats->PSRF), k);
7485 (*sumPSRF) += theStats->PSRF;
7486 (*numPSRFSamples)++;
7487 if (theStats->PSRF > *maxPSRF)
7488 (*maxPSRF) = theStats->PSRF;
7493 MrBayesPrint (" *");
7496 MrBayesPrintf (fpVstat, "\n");
7497 MrBayesPrint ("\n");
7501 /* PrintSumtTaxaInfo: Print information on pruned and absent taxa */
7502 void PrintSumtTaxaInfo (void)
7504 int i, j, lineWidth, numExcludedTaxa, len;
7507 /* print out information on absent taxa */
7508 numExcludedTaxa = 0;
7509 for (i=0; i<numTaxa; i++)
7510 if (sumtParams.absentTaxa[i] == YES)
7513 if (numExcludedTaxa > 0)
7515 if (numExcludedTaxa == 1)
7516 MrBayesPrint ("%s The following taxon was absent from trees:\n", spacer);
7518 MrBayesPrint ("%s The following %d taxa were absent from trees:\n", spacer, numExcludedTaxa);
7519 MrBayesPrint ("%s ", spacer);
7521 for (i=0; i<numTaxa; i++)
7523 if (sumtParams.absentTaxa[i] == YES)
7526 strcpy (tempStr, taxaNames[i]);
7527 len = (int) strlen(tempStr);
7531 MrBayesPrint ("\n%s ", spacer);
7534 if (numExcludedTaxa == 1)
7535 MrBayesPrint ("%s\n", tempStr);
7536 else if (numExcludedTaxa == 2 && j == 1)
7537 MrBayesPrint ("%s ", tempStr);
7538 else if (j == numExcludedTaxa)
7539 MrBayesPrint ("and %s\n", tempStr);
7541 MrBayesPrint ("%s, ", tempStr);
7544 MrBayesPrint ("\n");
7547 /* print out information on pruned taxa */
7548 numExcludedTaxa = 0;
7549 for (i=0; i<numTaxa; i++)
7550 if (taxaInfo[i].isDeleted == YES && sumtParams.absentTaxa[i] == NO)
7553 if (numExcludedTaxa > 0)
7555 if (numExcludedTaxa == 1)
7556 MrBayesPrint ("%s The following taxon was pruned from trees:\n", spacer);
7558 MrBayesPrint ("%s The following %d taxa were pruned from trees:\n", spacer, numExcludedTaxa);
7559 MrBayesPrint ("%s ", spacer);
7561 for (i=0; i<numTaxa; i++)
7563 if (taxaInfo[i].isDeleted == YES && sumtParams.absentTaxa[i] == NO)
7566 strcpy (tempStr, taxaNames[i]);
7567 len = (int) strlen(tempStr);
7571 MrBayesPrint ("\n%s ", spacer);
7574 if (numExcludedTaxa == 1)
7575 MrBayesPrint ("%s\n", tempStr);
7576 else if (numExcludedTaxa == 2 && j == 1)
7577 MrBayesPrint ("%s ", tempStr);
7578 else if (j == numExcludedTaxa)
7579 MrBayesPrint ("and %s\n", tempStr);
7581 MrBayesPrint ("%s, ", tempStr);
7584 MrBayesPrint ("\n");
7589 /* Range: Determine range for a vector of MrBFlt values */
7590 void Range (MrBFlt *vals, int nVals, MrBFlt *min, MrBFlt *max)
7592 SortMrBFlt (vals, 0, nVals-1);
7595 *max = vals[nVals-1];
7599 /* ResetTaxonSet: Reset included taxa and local outgroup number */
7600 void ResetTaxonSet (void)
7604 /* reset numLocalTaxa and localOutGroup */
7607 for (i=j=0; i<numTaxa; i++)
7609 if (taxaInfo[i].isDeleted == NO)
7611 if (i == outGroupNum)
7612 localOutGroup = numLocalTaxa;
7619 void ResetTranslateTable (void)
7623 for (i=0; i<numTranslates; i++)
7625 free (transFrom[i]);
7633 isTranslateDef = NO;
7634 isTranslateDiff = NO;
7638 int ShowConPhylogram (FILE *fp, PolyTree *t, int screenWidth)
7640 int i, j, k, nLines, from, to, treeWidth=0, barLength, printExponential,
7641 precision, width, newPos, curPos, nTimes, numSpaces, maxLabelLength;
7642 char *printLine, *markLine, temp[30], *label;
7643 MrBFlt scale, f, scaleBar;
7646 /* set max label length */
7647 maxLabelLength = 20;
7649 /* allocate space for label, printLine and markLine */
7650 printLine = (char *) SafeCalloc ((2*screenWidth+2),sizeof(char));
7651 label = (char *) SafeCalloc (maxLabelLength+1, sizeof(char));
7652 if (!printLine || !label)
7654 markLine = printLine + screenWidth + 1;
7656 /* calculate scale */
7659 for (i=t->nNodes-2; i>=0; i--)
7661 p = t->allDownPass[i];
7662 /* find distance to root in relevant units */
7663 if (sumtParams.isClock == YES && sumtParams.isCalibrated == NO)
7664 p->f = t->root->depth - p->depth;
7665 else if (sumtParams.isClock == YES && sumtParams.isCalibrated == YES)
7666 p->f = t->root->age - p->age;
7668 p->f = p->anc->f + p->length;
7669 if (p->left == NULL)
7671 f = p->f / (screenWidth - Label(p,YES,NULL,maxLabelLength) - 2);
7675 treeWidth = screenWidth - Label(p,YES,NULL,maxLabelLength) - 2;
7680 /* calculate x coordinates */
7681 for (i=0; i<t->nNodes; i++)
7683 p = t->allDownPass[i];
7684 p->x = (int) (0.5 + (p->f / scale));
7687 /* calculate y coordinates and lines to print */
7688 for (i=nLines=0; i<t->nNodes; i++)
7690 p = t->allDownPass[i];
7691 if (p->left != NULL)
7694 for (q=p->left->sib; q->sib!=NULL; q=q->sib)
7696 p->y = (int) (0.5 + ((p->left->y + q->y) / 2.0));
7706 /* print tree line by line */
7708 for (i=0; i<nLines; i++)
7710 MrBayesPrint ("%s ", spacer);
7711 /* empty printLine */
7712 for (j=0; j<screenWidth; j++)
7718 for (j=0; j<t->nNodes; j++)
7720 p = t->allDownPass[j];
7724 /* this branch should be printed */
7728 /* this is the root of the whole tree */
7729 printLine[p->x] = '+';
7733 /* this is an ordinary branch */
7736 for (k=from+1; k<=to; k++)
7738 if (p == p->anc->left)
7740 if (markLine[from] == 0)
7741 printLine[from] = '/';
7743 printLine[from] = '|';
7746 else if (p->sib == NULL)
7748 if (markLine[from] == 1)
7749 printLine[from] = '\\';
7751 printLine[from] = '|';
7757 printLine[to] = '+';
7759 printLine[to] = '|';
7763 /* add label if the branch is terminal */
7764 Label(p,YES,label,maxLabelLength);
7765 sprintf (printLine+to+2,"%s", label);
7770 /* check for cross branches */
7771 for (j=0; j<screenWidth; j++)
7773 if (markLine[j] >= 1 && printLine[j] == ' ')
7776 MrBayesPrintf (fp, "%s\n",printLine);
7780 k = (int) (floor (log10 (scale * 80)));
7781 scaleBar = pow (10, k);
7782 barLength = (int) (scaleBar / scale);
7788 else if (barLength > 40)
7793 else if (barLength > 16)
7799 if (t->isClock == YES)
7801 MrBayesPrint ("%s ", spacer);
7802 for (i=0; i<treeWidth; i++)
7804 nTimes = (int) (treeWidth / (scaleBar / scale));
7805 for (i=0; i<=nTimes; i++)
7806 printLine[treeWidth - (int)(i*scaleBar/scale)] = '|';
7807 MrBayesPrint ("%s\n", printLine);
7809 MrBayesPrint ("%s ", spacer);
7810 f = treeWidth * scale;
7811 if (f >= 1000.0 || f < 0.10)
7813 printExponential = YES;
7818 printExponential = NO;
7819 precision = 2 - (int) (log10 (f));
7823 f = nTimes * scaleBar;
7824 while (curPos < treeWidth)
7826 /* print the number */
7827 if (printExponential == YES)
7828 sprintf (temp, "%.2e", f);
7830 sprintf (temp, "%.*lf", precision, f);
7832 /* room to print ? if so, print */
7833 width = (int) strlen (temp);
7834 newPos = treeWidth - (int) (nTimes * (scaleBar / scale));
7835 numSpaces = newPos - width / 2 - curPos;
7836 if (numSpaces >= 0 || (numSpaces >= -2 && curPos == 0))
7838 while (numSpaces > 0)
7840 printLine[curPos++] = ' ';
7843 for (i=0; temp[i]!='\0'; i++)
7844 printLine[curPos++] = temp[i];
7847 /* get new number */
7852 MrBayesPrint ("%s\n", printLine);
7854 if (sumtParams.isCalibrated == YES)
7855 MrBayesPrint ("\n%s [User-defined time units]\n\n", spacer);
7857 MrBayesPrint ("\n%s [Expected changes per site]\n\n", spacer);
7861 MrBayesPrintf (fp, "%s |", spacer);
7862 for (i=0; i<barLength-1; i++)
7863 MrBayesPrintf (fp, "-");
7864 MrBayesPrintf (fp, "| %1.3lf expected changes per site\n\n", scaleBar);
7873 int ShowConTree (FILE *fp, PolyTree *t, int screenWidth, int showSupport)
7875 int i, j, k, treeWidth, minBranchLength, maxWidth, isTreeDivided,
7876 printWidth, nLines, nodesToBePrinted, from, to, maxLabelLength,
7878 char *printLine, *markLine, temp[20], *label;
7879 PolyNode *p=NULL, *q;
7881 maxLength = 20; /* max length of label */
7882 minBranchLength = 5; /* min length of branch in tree */
7885 /* allocate space for printLine, markLine and label */
7886 printLine = (char *) SafeCalloc (maxLength+1+(2*screenWidth+2),sizeof(char));
7889 markLine = printLine + screenWidth + 1;
7890 label = markLine + screenWidth + 1;
7892 /* get fresh internal node indices */
7893 k = t->nNodes - t->nIntNodes;
7894 for (i=0; i<t->nIntNodes; i++)
7896 p = t->intDownPass[i];
7900 /* calculate max length of labels including taxon index number */
7902 for (i=0; i<t->nNodes; i++)
7904 p = t->allDownPass[i];
7905 if (p->left == NULL)
7907 j = Label(p,YES,NULL,maxLength);
7908 if (j > maxLabelLength)
7913 /* make sure label can hold an interior node index number */
7914 j = (int) (3.0 + log10((MrBFlt)t->nNodes));
7915 maxLabelLength = (maxLabelLength > j ? maxLabelLength : j);
7917 /* calculate remaining screen width for tree
7918 and maxWidth in terms of branches */
7919 treeWidth = screenWidth - maxLabelLength - 1;
7920 maxWidth = treeWidth / minBranchLength;
7922 /* unmark whole tree */
7923 for (i=0; i<t->nNodes; i++)
7924 t->allDownPass[i]->mark = 0;
7925 nodesToBePrinted = t->nNodes;
7927 while (nodesToBePrinted > 0)
7929 /* count depth of nodes in unprinted tree */
7930 for (i=0; i<t->nNodes; i++)
7932 p = t->allDownPass[i];
7933 if (p->mark == 0) /* the node has not been printed yet */
7936 /* if it is an interior node in the tree that will be printed
7937 compute the depth of the node */
7938 if (p->left != NULL && p->left->mark == 0)
7940 for (q = p->left; q!=NULL; q=q->sib)
7946 /* break when root of print subtree has been found */
7947 if (p->x >= maxWidth)
7953 /* if internal node then find largest nonprinted subtree among descendant nodes */
7956 for (q=p->left; q!=NULL; q=q->sib)
7958 if (q->x == p->x - 1 && q->mark == 0)
7961 MrBayesPrintf (fp, "%s Subtree rooted at node %d:\n\n", spacer, p->index);
7962 isTreeDivided = YES;
7964 else if (isTreeDivided == YES)
7965 MrBayesPrintf (fp, "%s Root part of tree:\n\n", spacer);
7967 /* mark subtree for printing and
7968 translate x coordinates from depth to position */
7972 printWidth = p->x + 1;
7974 p->x = (int) (treeWidth - 0.5 - ((treeWidth - 1) * (p->x / (MrBFlt) printWidth)));
7975 for (i=t->nNodes-2; i>=0; i--)
7977 p = t->allDownPass[i];
7978 if (p->mark == 0 && p->anc->mark == 1)
7981 p->x = (int) (treeWidth - 0.5 - ((treeWidth - 1) * (p->x / (MrBFlt) printWidth)));
7985 /* calculate y coordinates of nodes to be printed and lines to print */
7986 for (i=nLines=0; i<t->nNodes; i++)
7988 p = t->allDownPass[i];
7991 if (p->left != NULL && p->left->mark == 1)
7994 for (q=p->left->sib; q->sib!=NULL; q=q->sib)
7996 p->y = (int) (0.5 + ((p->left->y + q->y) / 2.0));
8007 /* print subtree line by line */
8008 for (i=0; i<nLines; i++)
8010 MrBayesPrintf (fp, "%s ", spacer);
8011 /* empty printLine */
8012 for (j=0; j<screenWidth; j++)
8018 for (j=0; j<t->nNodes; j++)
8020 p = t->allDownPass[j];
8021 if (p->mark != 1 || p->y != i)
8024 /* this branch should be printed
8025 add label if the branch is terminal in tree to be printed */
8026 if (p->left == NULL)
8028 Label (p,YES,label,maxLength);
8029 sprintf (printLine+treeWidth+1,"%s", label);
8031 else if (p->left->mark == 2)
8032 sprintf (printLine+treeWidth+1,"(%d)", p->index);
8037 /* this is the root of the whole tree */
8038 printLine[p->x] = '+';
8041 else if (p->anc->mark == 0)
8043 /* this is a root of a subtree
8044 this branch will have to be printed again so do
8045 not decrease nodesToBePrinted */
8048 for (k=from; k<to; k++)
8050 printLine[to] = '+';
8051 if (showSupport == YES)
8052 sprintf (temp, "%d", (int) (p->support*100.0 + 0.5));
8055 from = (int)(from + 1.5 + ((to - from - 1 - strlen(temp)) / 2.0));
8056 for (k=0; temp[k]!='\0'; k++)
8057 printLine[from++] = temp[k];
8061 /* this is an ordinary branch */
8064 for (k=from+1; k<=to; k++)
8066 if (p == p->anc->left)
8068 printLine[from] = '/';
8071 else if (p->sib == NULL)
8073 printLine[from] = '\\';
8076 if (p->left!=NULL && p->left->mark!=2)
8078 printLine[to] = '+';
8079 if (showSupport == YES)
8080 sprintf (temp, "%d", (int) (p->support*100.0 + 0.5));
8083 from = (int)(from + 1.5 + ((to - from - 1 - strlen(temp)) / 2.0));
8084 for (k=0; temp[k]!='\0'; k++)
8085 printLine[from++] = temp[k];
8091 /* check for cross branches */
8092 for (j=0; j<treeWidth; j++)
8094 if (markLine[j] == 1 && printLine[j] == ' ')
8098 MrBayesPrintf (fp, "%s\n",printLine);
8101 /* mark printed branches */
8102 for (i=0; i<t->nNodes; i++)
8104 p = t->allDownPass[i];
8109 else if (p->anc->mark == 0)
8110 p->mark = 0; /* this branch will have to be printed again */
8116 } /* next subtree */
8124 void ShowParts (FILE *fp, BitsLong *p, int nTaxaToShow)
8127 BitsLong x, y, bitsLongOne;
8131 for (i=0; i<nTaxaToShow; i++)
8133 y = p[i / nBitsInALong];
8134 x = bitsLongOne << (i % nBitsInALong);
8136 MrBayesPrintf (fp, ".");
8138 MrBayesPrintf (fp, "*");
8143 void ShowSomeParts (FILE *fp, BitsLong *p, int offset, int nTaxaToShow)
8146 BitsLong x, y, bitsLongOne;
8150 for (i=offset; i<offset+nTaxaToShow; i++)
8152 y = p[i / nBitsInALong];
8153 x = bitsLongOne << (i % nBitsInALong);
8155 MrBayesPrintf (fp, ".");
8157 MrBayesPrintf (fp, "*");
8162 void SortPartCtr (PartCtr **item, int left, int right)
8165 PartCtr *tempPartCtr;
8169 assert (right >= 0);
8173 x = item[(left+right)/2]->totCount;
8176 while (item[i]->totCount > x && i < right)
8178 while (x > item[j]->totCount && j > left)
8182 tempPartCtr = item[i];
8184 item[j] = tempPartCtr;
8191 SortPartCtr (item, left, j);
8193 SortPartCtr (item, i, right);
8197 void SortTerminalPartCtr (PartCtr **item, int len)
8199 register int i, j, maxCount;
8202 maxCount = item[0]->totCount;
8204 /* put root first */
8205 for (i=0; item[i]->totCount == maxCount; i++)
8206 if (NumBits(item[i]->partition, sumtParams.BitsLongsNeeded) == sumtParams.numTaxa)
8216 /* then find terminals in index order */
8217 for (i=1; i<=sumtParams.numTaxa; i++)
8219 for (j=i; item[j]->totCount == maxCount && j<len; j++)
8220 if (NumBits(item[j]->partition, sumtParams.BitsLongsNeeded) == 1 &&
8221 FirstTaxonInPartition(item[j]->partition, sumtParams.BitsLongsNeeded) == i-1)
8234 void SortTreeCtr (TreeCtr **item, int left, int right)
8237 TreeCtr *tempTreeCtr;
8242 x = item[(left+right)/2]->count;
8245 while (item[i]->count > x && i < right)
8247 while (x > item[j]->count && j > left)
8251 tempTreeCtr = item[i];
8253 item[j] = tempTreeCtr;
8260 SortTreeCtr (item, left, j);
8262 SortTreeCtr (item, i, right);
8266 /* StoreSumtTree: Store tree in treeList in packed format */
8267 int StoreSumtTree (PackedTree *treeList, int index, PolyTree *t)
8269 int orderLen, numBrlens;
8271 assert (treeList[index].brlens == NULL);
8272 assert (treeList[index].order == NULL);
8274 /* get tree dimensions */
8275 numBrlens = t->nNodes - 1;
8276 orderLen = t->nIntNodes - 1;
8278 /* allocate space */
8279 treeList[index].brlens = (MrBFlt *) SafeCalloc (numBrlens, sizeof(MrBFlt));
8280 treeList[index].order = (int *) SafeCalloc (orderLen, sizeof(MrBFlt));
8281 if (!treeList[index].order || !treeList[index].brlens)
8283 MrBayesPrint ("%s Could not store packed representation of tree '%s'\n", spacer, t->name);
8288 if (t->isRooted == YES)
8289 StoreRPolyTree (t, treeList[index].order, treeList[index].brlens);
8291 StoreUPolyTree (t, treeList[index].order, treeList[index].brlens);
8297 /* TreeCtrUppass: extract TreeCtr nodes in uppass sequence */
8298 void TreeCtrUppass (TreeCtr *r, TreeCtr **uppass, int *index)
8302 uppass[(*index)++] = r;
8304 TreeCtrUppass (r->left, uppass, index);
8305 TreeCtrUppass (r->right, uppass, index);
8312 int i, num, nInSets[5];
8313 MrBFlt treeProb, cumTreeProb;
8317 /* check if we need to do this */
8318 if (sumtParams.calcTreeprobs == NO)
8321 MrBayesPrint ("%s Calculating tree probabilities...\n\n", spacer);
8323 /* allocate space for tree counters and trees */
8324 trees = (TreeCtr **) SafeCalloc ((size_t)numUniqueTreesFound, sizeof(TreeCtr *));
8325 theTree = AllocateTree (sumtParams.numTaxa);
8326 if (!trees || !theTree)
8328 MrBayesPrint ("%s Problem allocating trees or theTree in TreeProb\n", spacer);
8334 TreeCtrUppass (treeCtrRoot, trees, &i);
8337 SortTreeCtr (trees, 0, numUniqueTreesFound-1);
8339 /* set basic params in receiving tree */
8340 theTree->isRooted = sumtParams.isRooted;
8341 if (theTree->isRooted)
8343 theTree->nNodes = 2 * sumtParams.numTaxa;
8344 theTree->nIntNodes = sumtParams.numTaxa - 1;
8348 theTree->nNodes = 2 * sumtParams.numTaxa - 2;
8349 theTree->nIntNodes = sumtParams.numTaxa - 2;
8352 /* show tree data */
8354 nInSets[0] = nInSets[1] = nInSets[2] = nInSets[3] = nInSets[4] = 0;
8355 for (num=0; num<numUniqueTreesFound; num++) /* loop over all of the trees that were found */
8357 /* get probability of tree */
8358 treeProb = (MrBFlt)trees[num]->count / (MrBFlt)sumtParams.numTreesSampled;
8359 cumTreeProb += treeProb;
8360 if (cumTreeProb >= 0.0 && cumTreeProb < 0.5)
8362 else if (cumTreeProb >= 0.5 && cumTreeProb < 0.9)
8364 else if (cumTreeProb >= 0.9 && cumTreeProb < 0.95)
8366 else if (cumTreeProb >= 0.95 && cumTreeProb < 0.99)
8371 /* draw tree to stdout */
8372 if (theTree->isRooted == YES)
8373 RetrieveRTopology (theTree, trees[num]->order);
8375 RetrieveUTopology (theTree, trees[num]->order);
8376 if (sumtParams.showSumtTrees == YES)
8378 MrBayesPrint ("\n%s Tree %d (p = %1.3lf, P = %1.3lf):\n\n", spacer, num+1, treeProb, cumTreeProb);
8382 /* draw tree to file */
8385 MrBayesPrintf (fpTrees, "[This file contains the trees that were found during the MCMC\n");
8386 MrBayesPrintf (fpTrees, "search, sorted by posterior probability. \"p\" indicates the\n");
8387 MrBayesPrintf (fpTrees, "posterior probability of the tree whereas \"P\" indicates the\n");
8388 MrBayesPrintf (fpTrees, "cumulative posterior probability.]\n\n");
8389 MrBayesPrintf (fpTrees, "begin trees;\n");
8390 MrBayesPrintf (fpTrees, " translate\n");
8391 for (i=0; i<sumtParams.numTaxa; i++)
8393 if (i == sumtParams.numTaxa - 1)
8394 MrBayesPrintf (fpTrees, " %2d %s;\n", i+1, sumtParams.taxaNames[i]);
8396 MrBayesPrintf (fpTrees, " %2d %s,\n", i+1, sumtParams.taxaNames[i]);
8399 MrBayesPrintf (fpTrees, " tree tree_%d [p = %1.3lf, P = %1.3lf] = [&W %1.6lf] ", num+1, treeProb, cumTreeProb, treeProb);
8400 WriteTopologyToFile (fpTrees, theTree->root->left, theTree->isRooted);
8401 MrBayesPrintf (fpTrees, ";\n");
8402 if (num == numUniqueTreesFound - 1)
8403 MrBayesPrintf (fpTrees, "end;\n");
8406 /* print out general information on credible sets of trees */
8407 i = nInSets[0] + nInSets[1] + nInSets[2] + nInSets[3] + nInSets[4];
8408 MrBayesPrint ("%s Credible sets of trees (%d tree%s sampled):\n", spacer, i, i > 1 ? "s" : "");
8411 MrBayesPrint ("%s 50 %% credible set contains %d trees\n", spacer, i);
8414 MrBayesPrint ("%s 90 %% credible set contains %d trees\n", spacer, i);
8417 MrBayesPrint ("%s 95 %% credible set contains %d trees\n", spacer, i);
8419 MrBayesPrint ("%s 99 %% credible set contains %d tree%s\n\n", spacer, i, i > 1 ? "s" : "");
8428 void WriteConTree (PolyNode *p, FILE *fp, int showSupport)
8433 if (p->anc->left == p)
8436 for (q = p->left; q != NULL; q = q->sib)
8438 if (q->anc->left != q) /* Note that q->anc always exists (it is p) */
8440 WriteConTree (q, fp, showSupport);
8442 if (p->left == NULL)
8444 if (sumtParams.brlensDef == YES)
8446 if (sumtParams.isClock == NO)
8447 fprintf (fp, "%d:%s", p->index+1, MbPrintNum(p->length));
8449 fprintf (fp, "%d:%s", p->index+1, MbPrintNum(p->anc->depth - p->depth));
8452 fprintf (fp, "%d", p->index+1);
8455 if (p->sib == NULL && p->anc != NULL)
8457 if (p->anc->anc != NULL)
8459 if (sumtParams.brlensDef == YES && showSupport == NO)
8461 if (sumtParams.isClock == NO)
8462 fprintf (fp, "):%s", MbPrintNum(p->anc->length));
8464 fprintf (fp, "):%s", MbPrintNum(p->anc->anc->depth - p->anc->depth));
8466 else if (sumtParams.brlensDef == NO && showSupport == YES)
8467 fprintf (fp, ")%1.3lf", p->anc->support);
8468 else if (sumtParams.brlensDef == YES && showSupport == YES)
8470 if (sumtParams.isClock == NO)
8471 fprintf (fp, ")%1.3lf:%s", p->anc->support, MbPrintNum(p->anc->length));
8473 fprintf (fp, ")%1.3lf:%s", p->anc->support, MbPrintNum(p->anc->anc->depth - p->anc->depth));
8484 /* WriteFigTreeConTree: Include rich information for each node in a consensus tree */
8485 void WriteFigTreeConTree (PolyNode *p, FILE *fp, PartCtr **treeParts)
8489 if (p->left == NULL)
8491 fprintf (fp, "%d", p->index+1);
8492 if (sumtParams.isClock == NO)
8493 PrintFigTreeNodeInfo(fp, treeParts[p->partitionIndex], p->length);
8494 else if (sumtParams.isCalibrated == YES)
8495 PrintFigTreeNodeInfo(fp, treeParts[p->partitionIndex], p->anc->age - p->age);
8497 PrintFigTreeNodeInfo(fp, treeParts[p->partitionIndex], p->anc->depth - p->depth);
8502 for (q = p->left; q != NULL; q = q->sib)
8504 WriteFigTreeConTree (q, fp, treeParts);
8509 if (p->partitionIndex >= 0 && p->partitionIndex < numUniqueSplitsFound)
8513 if (sumtParams.isClock == YES)
8514 PrintFigTreeNodeInfo(fp,treeParts[p->partitionIndex], -1.0);
8516 else if (sumtParams.isClock == NO)
8517 PrintFigTreeNodeInfo(fp, treeParts[p->partitionIndex], p->length);
8518 else if (sumtParams.isCalibrated == YES)
8519 PrintFigTreeNodeInfo(fp, treeParts[p->partitionIndex], p->anc->age - p->age);
8521 PrintFigTreeNodeInfo(fp, treeParts[p->partitionIndex], p->anc->depth - p->depth);