/* * MrBayes 3.1.2 * * copyright 2002-2005 * * John P. Huelsenbeck * Section of Ecology, Behavior and Evolution * Division of Biological Sciences * University of California, San Diego * La Jolla, CA 92093-0116 * * johnh@biomail.ucsd.edu * * Fredrik Ronquist * School of Computational Science * Florida State University * Tallahassee, FL 32306-4120 * * ronquist@csit.fsu.edu * * This program is free software; you can redistribute it and/or * modify it under the terms of the GNU General Public License * as published by the Free Software Foundation; either version 2 * of the License, or (at your option) any later version. * * This program is distributed in the hope that it will be useful, * but WITHOUT ANY WARRANTY; without even the implied warranty of * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * GNU General Public License for more details (www.gnu.org). * */ #include #include #include #include #include #include #include "mb.h" #include "globals.h" #include "command.h" #include "bayes.h" #include "mbmath.h" #include "sumt.h" #include "mcmc.h" #if defined(__MWERKS__) #include "SIOUX.h" #endif typedef struct stnode { struct stnode *left, *right, *anc; int memoryIndex, index, upDateCl, upDateTi, marked, x, y, scalerNode, taxonName; long scalersSet, clSpace, tiSpace; char label[100]; MrBFlt length; } SumtNode; #define MAX_PARTITIONS 10000 #define MAX_TREES 1000 #undef DEBUG_CONTREE /* local prototypes */ int AddTreeToList (int whichList); int AllocBits (int n); void AssignIntPart (long *lft, long *rht, long *p); void AssignTipPart (int n, long *p); int BrlensVals (int treeNo, char *s, int longestLineLength, int lastTreeBlockBegin, int lastTreeBlockEnd); void CalculateTreeToTreeDistance (int *lst[2], MrBFlt *lngs[2], int nnds, int hasBrlens, MrBFlt *d1, MrBFlt *d2, MrBFlt *d3); int CheckSumtSpecies (void); int ConTree (void); int DerootSumtTree (SumtNode *p, int n, int outGrp); int FindParts (void); int FindTree (void); void FinishSumtTree (SumtNode *p, int *i, int isThisTreeRooted); int FirstTaxonInPartition (long *partition, int length); int FreeBits (void); void GetConDownPass (PolyNode **downPass, PolyNode *p, int *i); int GetPartitions (void); void GetSumtDownPass (SumtNode *p, SumtNode **dp, int *i); void GetSumtToken (int *tokenType); int OpenBrlensFile (int treeNo); int OpenComptFiles (void); int OpenSumtFiles (int treeNo); int PartFinder (long *p, MrBFlt bl, int *partID); int PrintBrlensToFile (void); void PrintParts (FILE *fp, long *p, int nTaxaToShow); int PruneSumt (void); int RealloateBits (void); int RealloateFullCompTrees (int whichList); int RealloateFullTrees (void); int ReorderParts (void); int RootSumtTree (SumtNode *p, int n, int out); void ShowBits (long *p, int nBitsToShow); void ShowConNodes (int nNodes, PolyNode *root); int ShowConTree (FILE *fp, int nNodes, PolyNode *root, int screenWidth, int showSupport); int ShowConPhylogram (FILE *fp, int nNodes, PolyNode *root, int screenWidth); void ShowParts (long *p, int nTaxaToShow); void ShowSumtNodes (SumtNode *p, int indent, int isThisTreeRooted); void SortIndParts (int *item, MrBFlt *assoc, int count, int descendingOrder); void SortIndParts2 (int *item, MrBFlt *assoc, int left, int right, int descendingOrder); int SortParts (int *item, int count); void SortParts2 (int *item, int left, int right); int SumtDex (SumtNode *p); int TreeProb (void); void WriteConTree (PolyNode *p, FILE *fp, int showSupport); void WriteTree (PolyNode *p, FILE *fp); extern int SafeFclose(FILE **); /* local (to this file) */ char *sumtTokenP, sumtToken[CMD_STRING_LENGTH]; int taxonLongsNeeded, numPartsAllocated, *numFoundOfThisPart, *numFoundOfThisPart1, *numFoundOfThisPart2, numTreesInLastBlock, **numFoundInRunOfPart, numTreePartsFound, numSumtTrees, numSumTreesSampled, numTranslates, whichTranslate, sumtBrlensDef, nextAvailableSumtNode, numSumtTaxa, isFirstSumtNode, foundSumtColon, *sumTaxaFound, numIncludedTaxa, isSumtTreeDefined, isSumtTreeRooted, *partOrigOrder, numAsterices, numTreeParts, *treePartNums, *fullTreePartIds, *numOfThisFullTree, numFullTreesAllocated, numFullTreesFound, *prunedTaxa, *absentTaxa, *firstPrunedTaxa, *firstAbsentTaxa, comparingFiles, fileNum, numCompTrees[2], numCompTreesSampled[2], numFullCompTreesFound[2], numFullCompTreesAllocated[2], *fullCompTreePartIds1, *fullCompTreePartIds2, numBrlens, printingBrlens, runIndex; long int *treeBits, *treePartsFound, *taxonMask; MrBFlt *aBrlens, *sBrlens, *treePartLengths, *fullCompTreePartLengths1, *fullCompTreePartLengths2, *brlens, *aWithinBrlens, *sWithinBrlens, *sumB, *sumsqB; SumtNode *pSumtPtr, *qSumtPtr, *sumtRoot, *sumtNodes; PolyNode *conNodes, *conRoot; FILE *fpParts, *fpCon, *fpTrees, *fpCompParts, *fpCompDists, *fpBrlens; int AddTreeToList (int whichList) { int i, *x; MrBFlt *y; if (numTreeParts == 0) { MrBayesPrint ("%s Too few tree partitions\n", spacer); return (ERROR); } if (numFullCompTreesFound[whichList] + 1 > numFullCompTreesAllocated[whichList]) { numFullCompTreesAllocated[whichList] += 500; if (RealloateFullCompTrees (whichList) == ERROR) return (ERROR); } if (whichList == 0) { x = &fullCompTreePartIds1[numFullCompTreesFound[0] * 2 * numTaxa]; for (i=0; i 1) { if (memAllocs[ALLOC_NUMINRUNOFPART] == YES) { MrBayesPrint ("%s numFoundOfThisPart not free in AllocBits\n", spacer); goto errorExit; } numFoundInRunOfPart = (int **) calloc ((size_t)(sumtParams.numRuns), sizeof (int *)); if (!numFoundInRunOfPart) { MrBayesPrint ("%s Problem allocating numFoundInRunOfPart (%d bytes)\n", spacer, sumtParams.numRuns * sizeof (int *)); goto errorExit; } for (i=0; i 1 && sumtBrlensDef == YES) { if (memAllocs[ALLOC_A_WITHIN_BRLENS] == YES) { MrBayesPrint ("%s aWithinBrlens not free in AllocBits\n", spacer); goto errorExit; } aWithinBrlens = (MrBFlt *)SafeMalloc((size_t) (MAX_PARTITIONS * sizeof(MrBFlt))); if (!aWithinBrlens) { MrBayesPrint ("%s Problem allocating aWithinBrlens (%d)\n", spacer, MAX_PARTITIONS * sizeof(MrBFlt)); goto errorExit; } memAllocs[ALLOC_A_WITHIN_BRLENS] = YES; if (memAllocs[ALLOC_S_WITHIN_BRLENS] == YES) { MrBayesPrint ("%s sWithinBrlens not free in AllocBits\n", spacer); goto errorExit; } sWithinBrlens = (MrBFlt *)SafeMalloc((size_t) (MAX_PARTITIONS * sizeof(MrBFlt))); if (!sWithinBrlens) { MrBayesPrint ("%s Problem allocating sWithinBrlens (%d)\n", spacer, MAX_PARTITIONS * sizeof(MrBFlt)); goto errorExit; } memAllocs[ALLOC_S_WITHIN_BRLENS] = YES; if (memAllocs[ALLOC_SUMB] == YES) { MrBayesPrint ("%s sumB not free in AllocBits\n", spacer); goto errorExit; } sumB = (MrBFlt *)SafeMalloc((size_t) (MAX_PARTITIONS * sizeof(MrBFlt))); if (!sumB) { MrBayesPrint ("%s Problem allocating sumB (%d bytes)\n", spacer, MAX_PARTITIONS * sizeof(MrBFlt)); goto errorExit; } memAllocs[ALLOC_SUMB] = YES; if (memAllocs[ALLOC_SUMSQB] == YES) { MrBayesPrint ("%s sumsqB not free in AllocBits\n", spacer); goto errorExit; } sumsqB = (MrBFlt *)SafeMalloc((size_t) (MAX_PARTITIONS * sizeof(MrBFlt))); if (!sumsqB) { MrBayesPrint ("%s Problem allocating sumsqB (%d bytes)\n", spacer, MAX_PARTITIONS * sizeof(MrBFlt)); goto errorExit; } memAllocs[ALLOC_SUMSQB] = YES; } if (memAllocs[ALLOC_TAXAFOUND] == YES) { MrBayesPrint ("%s sumTaxaFound not free in AllocBits\n", spacer); goto errorExit; } sumTaxaFound = NULL; sumTaxaFound = (int *)SafeMalloc((size_t) (numTaxa * sizeof(int))); if (!sumTaxaFound) { MrBayesPrint ("%s Problem allocating sumTaxaFound (%d)\n", spacer, numTaxa * sizeof(int)); goto errorExit; } for (i=0; i 1) i = 4 * numTaxa; else i = 2 * numTaxa; prunedTaxa = (int *)SafeMalloc((size_t) (i * sizeof(int))); if (!prunedTaxa) { MrBayesPrint ("%s Problem allocating prunedTaxa (%d)\n", spacer, i * sizeof(long)); goto errorExit; } absentTaxa = prunedTaxa + numTaxa; if (sumtParams.numRuns > 1) { firstPrunedTaxa = prunedTaxa + 2 * numTaxa; firstAbsentTaxa = prunedTaxa + 3 * numTaxa; } memAllocs[ALLOC_PRUNEINFO] = YES; if (memAllocs[ALLOC_SUMTTREE] == YES) { MrBayesPrint ("%s sumtNodes not free in AllocBits\n", spacer); goto errorExit; } sumtNodes = NULL; sumtNodes = (SumtNode *)SafeMalloc((size_t) (2 * numTaxa * sizeof(SumtNode))); if (!sumtNodes) { MrBayesPrint ("%s Problem allocating sumtNodes (%d)\n", spacer, 2 * numTaxa * sizeof(SumtNode)); goto errorExit; } for (i=0; i<2*numTaxa; i++) { sumtNodes[i].left = sumtNodes[i].right = sumtNodes[i].anc = NULL; sumtNodes[i].memoryIndex = i; sumtNodes[i].length = 0.0; sumtNodes[i].marked = NO; sumtNodes[i].index = 0; } sumtRoot = NULL; isSumtTreeDefined = NO; memAllocs[ALLOC_SUMTTREE] = YES; /* initialize */ for (i=0; i<2*n*taxonLongsNeeded; i++) treeBits[i] = 0; for (i=0; i 1) { for (i=0; i 1) { for (i=0; i longIntegerSize*8*(offSet+1)) offSet++; x = 1 << (i - offSet*longIntegerSize*8); y = taxonMask[offSet]; taxonMask[offSet] = x | y; } for (i=0; i 1) { for (i=0; i= sumtParams.brlensFreqDisplay; i++) ; numBrlens = i; if (numBrlens < 1) { MrBayesPrint ("%s No branch lengths above the frequency to display (%lf).", spacer, sumtParams.brlensFreqDisplay); return ERROR; } /* Open brlens file */ if (OpenBrlensFile(treeNo) == ERROR) return ERROR; /* allocate space for brlens */ brlens = (MrBFlt *) SafeMalloc (numBrlens * sizeof(MrBFlt)); if (brlens == NULL) { SafeFclose (&fpBrlens); return ERROR; } for (i=0; i 1 && sumtParams.numTrees == 1) sprintf (fileName, "%s.run%d.t", sumtParams.sumtFileName, runNo+1); else if (sumtParams.numRuns == 1 && sumtParams.numTrees > 1) sprintf (fileName, "%s.tree%d.t", sumtParams.sumtFileName, treeNo+1); else if (sumtParams.numRuns > 1 && sumtParams.numTrees > 1) sprintf (fileName, "%s.tree%d.run%d.t", sumtParams.sumtFileName, treeNo+1, runNo+1); /* open binary file */ if ((fp = OpenBinaryFileR (fileName)) == NULL) { SafeFclose (&fpBrlens); free (brlens); return ERROR; } /* find length of longest line */ longestLineLength = LongestLine (fp); longestLineLength += 10; /* allocate a string long enough to hold a line */ if (runNo == 0) s = (char *) SafeMalloc (sizeof (char) * longestLineLength); else { free (s); s = (char *) SafeMalloc (sizeof (char) * longestLineLength); } if (!s) { free (brlens); SafeFclose (&fpBrlens); return ERROR; } /* close binary file */ SafeFclose (&fp); /* open text file */ if ((fp = OpenTextFileR (fileName)) == NULL) { SafeFclose (&fpBrlens); free (brlens); free (s); return ERROR; } /* Check file for appropriate blocks. We want to find the last tree block in the file and start from there. */ foundBegin = inTreeBlock = inSumtComment = NO; lineNum = lastTreeBlockBegin = lastTreeBlockEnd = 0; while (fgets (s, longestLineLength, fp) != NULL) { sumtTokenP = &s[0]; do { GetSumtToken (&tokenType); if (IsSame("[", sumtToken) == SAME) inSumtComment = YES; if (IsSame("]", sumtToken) == SAME) inSumtComment = NO; if (inSumtComment == NO) { if (foundBegin == YES) { if (IsSame("Trees", sumtToken) == SAME) { inTreeBlock = YES; foundBegin = NO; lastTreeBlockBegin = lineNum; } } else { if (IsSame("Begin", sumtToken) == SAME) { foundBegin = YES; } else if (IsSame("End", sumtToken) == SAME) { if (inTreeBlock == YES) { inTreeBlock = NO; lastTreeBlockEnd = lineNum; } } } } } while (*sumtToken); lineNum++; } /* Now fast rewind tree file */ (void)fseek(fp, 0L, 0); /* ...and fast forward to beginning of last tree block. */ for (i=0; ileft == NULL && q->right == NULL) || (q->anc == NULL && isSumtTreeRooted == NO)) { if (CheckString (q->label, taxaNames, &whichTaxon) == ERROR) { MrBayesPrint ("%s Could not find taxon %s in original list of taxa\n", spacer, q->label); goto errorExit; } whichTaxon--; sumTaxaFound[whichTaxon] = YES; } } /* now, print out which taxa are not in list */ numNotFound = 0; for (i=0; ix counts number of subtended terminals make sure conRoot->left is in outgroup */ conRoot = &conNodes[numIncludedTaxa]; conRoot->anc = conRoot->sib = NULL; conRoot->x = numIncludedTaxa; j = FirstTaxonInPartition (outgroupPartition, taxonLongsNeeded); conRoot->left = cp = &conNodes[j]; cp->anc = conRoot; cp->x = 1; for (i=0; isib = &conNodes[i]; cp = cp->sib; cp->anc = conRoot; cp->x = 1; } } cp->sib = NULL; /* Resolve bush according to partitions. Partitions may include incompatible ones. Partitions should be sorted from most frequent to least frequent for quit test to work when a 50% majority rule tree is requested. */ nextConNode = numIncludedTaxa + 1; if (isSumtTreeRooted == YES) targetNode = 2 * numIncludedTaxa - 2; else targetNode = 2 * numIncludedTaxa - 3; numTerminalsEncountered = 0; for (i=0; i targetNode && numTerminalsEncountered == numIncludedTaxa) break; freq = (MrBFlt)numFoundOfThisPart[i]/ (MrBFlt)(sumtParams.numRuns * numSumTreesSampled); if (freq < 0.50 && !strcmp(sumtParams.sumtConType, "Halfcompat")) break; /* get partition */ partition = &treePartsFound[i*taxonLongsNeeded]; /* flip bits if necessary */ /* This code is needed if single outgroup is indexed incorrectly or if the partition defines a clade in a multispecies outgroup but the indexing is reversed. Note that bits should not be flipped for rooted trees. */ if (isSumtTreeRooted == NO) { if (!IsPartNested(partition, ingroupPartition, taxonLongsNeeded) && !IsPartCompatible(partition, ingroupPartition, taxonLongsNeeded)) FlipBits(partition, taxonLongsNeeded, mask); } /* count bits in this partition */ for (j=nBits=0; j 1) { /* find anc of partition */ j = FirstTaxonInPartition (partition, taxonLongsNeeded); for (cp = &conNodes[j]; cp!=NULL; cp = cp->anc) if (cp->x > nBits) break; /* do not include if incompatible with ancestor or any of descendants do not check terminals or root because it is redundant and partitions have not been set for those */ isCompat = YES; if (cp->anc != NULL && !IsPartCompatible(partition, cp->partition, taxonLongsNeeded)) isCompat = NO; for (q=cp->left; q!=NULL; q=q->sib) { if (q->x > 1 && !IsPartCompatible(q->partition, partition, taxonLongsNeeded)) isCompat = NO; if (isCompat == NO) break; } if (isCompat == NO) continue; /* set new node */ q = &conNodes[nextConNode++]; if (sumtBrlensDef == YES) q->length = aBrlens[i]; else q->length = 0.0; q->support = freq * 100; q->x = nBits; q->partition = partition; /* go through descendants of anc */ ql = pl = NULL; for (r=cp->left; r!=NULL; r=r ->sib) { /* test if r is in the new partition or not */ if ((r->x > 1 && IsPartNested(r->partition, partition, taxonLongsNeeded)) || (r->x == 1 && (partition[r->index / nBitsInALong] & (1 << (r->index % nBitsInALong))) != 0)) { /* r is in the partition */ if (ql == NULL) q->left = r; else ql->sib = r; ql = r; r->anc = q; } else { /* r is not in the partition */ if (pl == NULL) cp->left = r; else pl->sib = r; pl = r; } } /* terminate new sib-node chain */ ql->sib = NULL; /* new node is last in old sib-node chain */ pl->sib = q; q->sib = NULL; q->anc = cp; } else /* singleton partition */ { j = FirstTaxonInPartition(partition, taxonLongsNeeded); q = &conNodes[j]; if (sumtBrlensDef == YES) q->length = aBrlens[i]; else q->length = 0.0; q->support = freq * 100; numTerminalsEncountered++; } } if (sumtParams.orderTaxa == YES) { /* rearrange tree so that terminals are in order */ /* first allocate space for downPass */ downPass = (PolyNode **) calloc (nextConNode, sizeof (PolyNode *)); if (!downPass) return ERROR; i = 0; GetConDownPass (downPass, conRoot, &i); /* label by minimum index */ for (i=0; ileft == NULL) { if (cp->index == localOutgroupNum) cp->x = -1; else cp->x = cp->index; } else { j = nextConNode; for (q=cp->left; q!=NULL; q=q->sib) { if (q->x < j) j = q->x; } cp->x = j; } } /* and rearrange */ for (i=0; ileft == NULL || cp->anc == NULL) continue; for (ql=NULL, q=cp->left; q->sib!=NULL; ql=q, q=q->sib) { for (rl=q, r=q->sib; r!=NULL; rl=r, r=r->sib) { if (r->x < q->x) { if (ql == NULL) cp->left = r; if (r == q->sib) /* swap adjacent q and r */ { if (ql != NULL) ql->sib = r; pl = r->sib; r->sib = q; q->sib = pl; } else /* swap separated q and r */ { if (ql != NULL) ql->sib = r; pl = r->sib; r->sib = q->sib; rl->sib = q; q->sib = pl; } pl = q; q = r; r = pl; } } } } } /* draw tree to stdout and fp */ MrBayesPrint ("\n%s Clade credibility values:\n\n", spacer); ShowConTree (stdout, nextConNode, conRoot, 80, YES); if (logToFile == YES) ShowConTree (logFileFp, nextConNode, conRoot, 80, YES); if (sumtBrlensDef == YES) { MrBayesPrint ("\n"); MrBayesPrint ("%s Phylogram:\n\n", spacer); ShowConPhylogram (stdout, nextConNode, conRoot, 80); if (logToFile == YES) ShowConPhylogram (logFileFp, nextConNode, conRoot, 80); } /* remove terminal taxon numbers */ for (i=0; i 0; j--) if (conNodes[i].label[j] == '(') break; conNodes[i].label[j-1] = '\0'; } MrBayesPrintf (fpCon, "begin trees;\n"); MrBayesPrintf (fpCon, " [Note: This tree contains information on the topology, \n"); MrBayesPrintf (fpCon, " branch lengths (if present), and the probability\n"); MrBayesPrintf (fpCon, " of the partition indicated by the branch.]\n"); if (!strcmp(sumtParams.sumtConType, "Halfcompat")) MrBayesPrintf (fpCon, " tree con_50_majrule = "); else MrBayesPrintf (fpCon, " tree con_all_compat = "); WriteConTree (conRoot, fpCon, YES); MrBayesPrintf (fpCon, ";\n"); if (sumtBrlensDef == YES) { MrBayesPrintf (fpCon, "\n"); MrBayesPrintf (fpCon, " [Note: This tree contains information only on the topology\n"); MrBayesPrintf (fpCon, " and branch lengths (mean of the posterior probability density).]\n"); if (!strcmp(sumtParams.sumtConType, "Halfcompat")) MrBayesPrintf (fpCon, " tree con_50_majrule = "); else MrBayesPrintf (fpCon, " tree con_all_compat = "); WriteConTree (conRoot, fpCon, NO); MrBayesPrintf (fpCon, ";\n"); } MrBayesPrintf (fpCon, "end;\n"); /* free memory and file pointers */ if (memAllocs[ALLOC_CONNODES] == YES) { free (conNodes); memAllocs[ALLOC_CONNODES] = NO; } if (memAllocs[ALLOC_OUTPART] == YES) { free (outgroupPartition); memAllocs[ALLOC_OUTPART] = NO; } if (downPass) free (downPass); return (NO_ERROR); errorExit: if (memAllocs[ALLOC_CONNODES] == YES) { free (conNodes); memAllocs[ALLOC_CONNODES] = NO; } if (memAllocs[ALLOC_OUTPART] == YES) { free (outgroupPartition); memAllocs[ALLOC_OUTPART] = NO; } if (downPass) free (downPass); return (ERROR); } int DerootSumtTree (SumtNode *p, int n, int outGrp) { int i, nNodes, isMarked, *localTaxaFound=NULL, sumtOut; MrBFlt tempBrLen; SumtNode **downPass=NULL, *lft, *m1, *m2, *um1, *um2, *out; # if 0 ShowTree (sumtRoot, YES, n); # endif /* check that we are not already unrooted */ if (isSumtTreeRooted == NO) { MrBayesPrint ("%s Tree is already unrooted\n", spacer); goto errorExit; } /* Find the outgroup number (0, 1, ..., n-1) and the number of nodes on the tree. The number of nodes may change later, if the user deleted taxa. */ sumtOut = outGrp; nNodes = 2 * n; /* allocate space for derooting tree */ downPass = (SumtNode **)SafeMalloc((size_t) (2 * numTaxa * sizeof(SumtNode *))); if (!downPass) { MrBayesPrint ("%s Could not allocate downPass\n", spacer); goto errorExit; } localTaxaFound = (int *)SafeMalloc((size_t) (numTaxa * sizeof(int))); if (!localTaxaFound) { MrBayesPrint ("%s Could not allocate localTaxaFound\n", spacer); goto errorExit; } for (i=0; imarked = NO; if (p->left == NULL && p->right == NULL && p->anc != NULL) { localTaxaFound[p->index] = YES; if (p->index == sumtOut) { p->marked = YES; isMarked = YES; } } } /* If we don't find the outgroup, it was deleted by the user. Instead, we will use the first undeleted taxon in the matrix that is found in the tree. */ if (isMarked == NO) { /* here we find the first undeleted taxon and designate it as the outgroup */ for (i=0; imarked = NO; if (p->left == NULL && p->right == NULL && p->anc != NULL) { if (p->index == sumtOut) { p->marked = YES; isMarked = YES; } } } /* if we still have not marked an outgroup, we have trouble */ if (isMarked == NO) { MrBayesPrint ("%s Could not find outgroup taxon\n", spacer); goto errorExit; } } /* now we mark the path from the outgroup tip to the root */ for (i=0; ileft != NULL && p->right != NULL) if (p->left->marked == YES || p->right->marked == YES) p->marked = YES; } /* now we rotate the tree until the outgroup is to the left or to the right of the root */ lft = sumtRoot->left; while (lft->left->index != sumtOut && lft->right->index != sumtOut) { if (lft->left->marked == YES && lft->right->marked == NO) { m1 = lft->left; um1 = lft->right; if (m1->left != NULL && m1->right != NULL) { if (m1->left->marked == YES) { m2 = m1->left; um2 = m1->right; } else { m2 = m1->right; um2 = m1->left; } lft->left = m2; lft->right = m1; m2->anc = m1->anc = lft; m1->left = um2; m1->right = um1; um1->anc = um2->anc = m1; m1->marked = NO; um1->length += m1->length; m2->length *= 0.5; m1->length = m2->length; } else { MrBayesPrint ("%s Rooting routine is lost (1)\n", spacer); goto errorExit; } } else if (lft->left->marked == NO && lft->right->marked == YES) { m1 = lft->right; um1 = lft->left; if (m1->left != NULL && m1->right != NULL) { if (m1->left->marked == YES) { m2 = m1->left; um2 = m1->right; } else { m2 = m1->right; um2 = m1->left; } lft->left = m1; lft->right = m2; m2->anc = m1->anc = lft; m1->left = um1; m1->right = um2; um1->anc = um2->anc = m1; m1->marked = NO; um1->length += m1->length; m2->length *= 0.5; m1->length = m2->length; } else { MrBayesPrint ("%s Rooting routine is lost (2)\n", spacer); goto errorExit; } } else { MrBayesPrint ("%s Rooting routine is lost (3)\n", spacer); goto errorExit; } } /* make certain outgroup is to the right of the root */ if (sumtRoot->left->left->index == sumtOut) { m1 = sumtRoot->left->left; m2 = sumtRoot->left->right; lft = sumtRoot->left; lft->left = m2; lft->right = m1; } /* now, take outgroup and make it point down */ m1 = sumtRoot; m2 = sumtRoot->left; lft = sumtRoot->left->left; out = sumtRoot->left->right; tempBrLen = out->length; if (tempBrLen < 0.0) tempBrLen = 0.0; lft->anc = out; out->left = lft; out->right = out->anc = NULL; lft->length += tempBrLen; m1->left = m1->right = m1->anc = NULL; m2->left = m2->right = m2->anc = NULL; sumtRoot = out; /* the tree is now unrooted */ isSumtTreeRooted = NO; /* reindex internal nodes of tree */ i = numTaxa; FinishSumtTree (sumtRoot, &i, NO); # if 0 ShowNodes (sumtRoot, 3, isSumtTreeRooted); ShowTree (sumtRoot, isSumtTreeRooted, n); # endif /* free memory */ free (downPass); free (localTaxaFound); return (NO_ERROR); errorExit: if (downPass) free (downPass); if (localTaxaFound) free (localTaxaFound); return(ERROR); } int DoCompareTree (void) { int i, j, k, n, lineTerm, longestLineLength, tokenType, foundBegin, inTreeBlock, lineNum, numTreesInBlock, numTreeBlocks, lastTreeBlockBegin, lastTreeBlockEnd, lineWidth, blockErrors, inSumtComment, numExcludedTaxa, longestLineLength1, longestLineLength2, xaxis, yaxis, starHolder[80], *treePartList[2], minNumTrees, screenWidth, screenHeigth, numY[60], nSamples, nCompPartitions, numIncludedTaxa1, numIncludedTaxa2, oldSumtBrlensDef, brlensDef; long int temporarySeed, len; long *x; MrBFlt xProb, yProb, xInc, yInc, xUpper, xLower, yUpper, yLower, *treeLengthList[2], *dT1=NULL, *dT2=NULL, *dT3=NULL, d1, d2, d3, meanY[60], xVal, yVal, minX, minY, maxX, maxY, sums[3]; char *s=NULL, tempStr[100], tempName[100], prCh; FILE *fp[2]; time_t curTime; # if defined (MPI_ENABLED) if (proc_id == 0) { # endif /* Are we comparing two files, or using sumt? */ comparingFiles = YES; fileNum = 0; oldSumtBrlensDef = sumtBrlensDef; /* we do not want to print brlens to file */ printingBrlens = NO; /* set file pointers to NULL */ fp[0] = fp[1] = NULL; /* Check that a data set has been read in. We check taxon names against those read in. */ if (defMatrix == NO) { MrBayesPrint ("%s A matrix must be specified before comparetree can be used\n", spacer); goto errorExit; } /* open binary file */ if ((fp[0] = OpenBinaryFileR(comptreeParams.comptFileName1)) == NULL) goto errorExit; if ((fp[1] = OpenBinaryFileR(comptreeParams.comptFileName2)) == NULL) goto errorExit; /* find out what type of line termination is used for file 1 */ lineTerm = LineTermType (fp[0]); if (lineTerm == LINETERM_MAC) MrBayesPrint ("%s Macintosh line termination for file 1\n", spacer); else if (lineTerm == LINETERM_DOS) MrBayesPrint ("%s DOS line termination for file 1\n", spacer); else if (lineTerm == LINETERM_UNIX) MrBayesPrint ("%s UNIX line termination for file 1\n", spacer); else { MrBayesPrint ("%s Unknown line termination for file 1\n", spacer); goto errorExit; } /* find out what type of line termination is used for file 2 */ lineTerm = LineTermType (fp[1]); if (lineTerm == LINETERM_MAC) MrBayesPrint ("%s Macintosh line termination for file 2\n", spacer); else if (lineTerm == LINETERM_DOS) MrBayesPrint ("%s DOS line termination for file 2\n", spacer); else if (lineTerm == LINETERM_UNIX) MrBayesPrint ("%s UNIX line termination for file 2\n", spacer); else { MrBayesPrint ("%s Unknown line termination for file 2\n", spacer); goto errorExit; } /* find length of longest line in either file */ longestLineLength1 = LongestLine (fp[0]); longestLineLength2 = LongestLine (fp[1]); if (longestLineLength1 > longestLineLength2) longestLineLength = longestLineLength1; else longestLineLength = longestLineLength2; MrBayesPrint ("%s Longest line length = %d\n", spacer, longestLineLength); longestLineLength += 10; /* allocate a string long enough to hold a line */ if (memAllocs[ALLOC_SUMTSTRING] == YES) { MrBayesPrint ("%s Comparetree string is already allocated\n", spacer); goto errorExit; } s = (char *)SafeMalloc((size_t) (longestLineLength * sizeof(char))); if (!s) { MrBayesPrint ("%s Problem allocating string for reading comparetree file\n", spacer); goto errorExit; } memAllocs[ALLOC_SUMTSTRING] = YES; /* close binary file */ SafeFclose (&fp[0]); SafeFclose (&fp[1]); /* read in data file 1 ***************************************************************************/ /* tell user we are ready to go */ MrBayesPrint ("%s Summarizing trees in file %s\n", spacer, comptreeParams.comptFileName1); /* open text file */ if ((fp[0] = OpenTextFileR(comptreeParams.comptFileName1)) == NULL) goto errorExit; /* Check file for appropriate blocks. We want to find the last tree block in the file and start from there. */ foundBegin = inTreeBlock = blockErrors = inSumtComment = NO; lineNum = numTreesInBlock = lastTreeBlockBegin = lastTreeBlockEnd = numTreeBlocks = numTreesInLastBlock = 0; while (fgets (s, longestLineLength, fp[0]) != NULL) { sumtTokenP = &s[0]; do { GetSumtToken (&tokenType); if (IsSame("[", sumtToken) == SAME) inSumtComment = YES; if (IsSame("]", sumtToken) == SAME) inSumtComment = NO; if (inSumtComment == NO) { if (foundBegin == YES) { if (IsSame("Trees", sumtToken) == SAME) { numTreesInBlock = 0; inTreeBlock = YES; foundBegin = NO; lastTreeBlockBegin = lineNum; } } else { if (IsSame("Begin", sumtToken) == SAME) { if (foundBegin == YES) { MrBayesPrint ("%s Found inappropriate \"Begin\" statement in file\n", spacer); blockErrors = YES; } foundBegin = YES; } else if (IsSame("End", sumtToken) == SAME) { if (inTreeBlock == YES) { numTreeBlocks++; inTreeBlock = NO; lastTreeBlockEnd = lineNum; } else { MrBayesPrint ("%s Found inappropriate \"End\" statement in file\n", spacer); blockErrors = YES; } numTreesInLastBlock = numTreesInBlock; } else if (IsSame("Tree", sumtToken) == SAME) { if (inTreeBlock == YES) { brlensDef = NO; for (j=0; s[j]!='\0'; j++) { if (s[j] == ':') { brlensDef = YES; break; } else if (s[j] == ',') break; } sumtBrlensDef = brlensDef; numTreesInBlock++; } else { MrBayesPrint ("%s Found a \"Tree\" statement that is not in a tree block\n", spacer); blockErrors = YES; } } } } } while (*sumtToken); lineNum++; } /* Now, check some aspects of the tree file, such as the number of tree blocks and whether they are properly terminated. */ if (inTreeBlock == YES) { MrBayesPrint ("%s Unterminated tree block in file %s. You probably need to\n", spacer, comptreeParams.comptFileName1); MrBayesPrint ("%s add a new line to the end of the file with \"End;\" on it.\n", spacer); goto errorExit; } if (inSumtComment == YES) { MrBayesPrint ("%s Unterminated comment in file %s\n", spacer, comptreeParams.comptFileName1); goto errorExit; } if (blockErrors == YES) { MrBayesPrint ("%s Found formatting errors in file %s\n", spacer, comptreeParams.comptFileName1); goto errorExit; } if (lastTreeBlockEnd < lastTreeBlockBegin) { MrBayesPrint ("%s Problem reading tree file %s\n", spacer, comptreeParams.comptFileName1); goto errorExit; } if (numTreesInLastBlock <= 0) { MrBayesPrint ("%s No trees were found in last tree block of file %s\n", spacer, comptreeParams.comptFileName1); goto errorExit; } if (sumtParams.sumtBurnIn > numTreesInLastBlock) { MrBayesPrint ("%s No trees are sampled as the burnin exceeds the number of trees in last block\n", spacer); MrBayesPrint ("%s Try setting burnin to a number less than %d\n", spacer, numTreesInLastBlock); goto errorExit; } /* tell the user that everything is fine */ if (numTreeBlocks == 1) MrBayesPrint ("%s Found one tree block in file \"%s\" with %d trees in last block\n", spacer, comptreeParams.comptFileName1, numTreesInLastBlock); else { MrBayesPrint ("%s Found %d tree blocks in file \"%s\" with %d trees in last block\n", spacer, numTreeBlocks, comptreeParams.comptFileName1, numTreesInLastBlock); MrBayesPrint ("%s Only the %d trees in last tree block will be summarized\n", spacer, numTreesInLastBlock); } /* Now we read the file for real. First, rewind file pointer to beginning of file... */ (void)fseek(fp[0], 0L, 0); /* ...and fast forward to beginning of last tree block. */ for (i=0; i 0) { if (numExcludedTaxa == 1) MrBayesPrint ("%s The following species was absent from trees:\n", spacer); else MrBayesPrint ("%s The following %d species were absent from trees:\n", spacer, numExcludedTaxa); MrBayesPrint ("%s ", spacer); j = lineWidth = 0; for (i=0; i 60) { MrBayesPrint ("\n%s ", spacer); lineWidth = 0; } if (numExcludedTaxa == 1) MrBayesPrint ("%s\n", tempStr); else if (numExcludedTaxa == 2 && j == 1) MrBayesPrint ("%s ", tempStr); else if (j == numExcludedTaxa) MrBayesPrint ("and %s\n", tempStr); else MrBayesPrint ("%s, ", tempStr); } } MrBayesPrint ("\n"); } /* print out information on pruned taxa */ numExcludedTaxa = 0; for (i=0; i 0) { if (numExcludedTaxa == 1) MrBayesPrint ("%s The following species was pruned from trees:\n", spacer); else MrBayesPrint ("%s The following %d species were pruned from trees:\n", spacer, numExcludedTaxa); MrBayesPrint ("%s ", spacer); j = lineWidth = 0; for (i=0; i 60) { MrBayesPrint ("\n%s ", spacer); lineWidth = 0; } if (numExcludedTaxa == 1) MrBayesPrint ("%s\n", tempStr); else if (numExcludedTaxa == 2 && j == 1) MrBayesPrint ("%s ", tempStr); else if (j == numExcludedTaxa) MrBayesPrint ("and %s\n", tempStr); else MrBayesPrint ("%s, ", tempStr); } } MrBayesPrint ("\n"); } /* tell user how many trees were successfully read */ MrBayesPrint ("%s Read %d trees from last tree block (sampling %d of them)\n", spacer, numSumtTrees, numSumTreesSampled); /* Check that at least one tree was read in. */ if (numSumTreesSampled <= 0) { MrBayesPrint ("%s No trees read in\n", spacer); goto errorExit; } /* read in data file 2 ***************************************************************************/ fileNum = 1; /* tell user we are ready to go */ MrBayesPrint ("%s Summarizing trees in file %s\n", spacer, comptreeParams.comptFileName2); /* open text file */ if ((fp[1] = OpenTextFileR(comptreeParams.comptFileName2)) == NULL) goto errorExit; /* Check file for appropriate blocks. We want to find the last tree block in the file and start from there. */ foundBegin = inTreeBlock = blockErrors = inSumtComment = NO; lineNum = numTreesInBlock = lastTreeBlockBegin = lastTreeBlockEnd = numTreeBlocks = numTreesInLastBlock = 0; while (fgets (s, longestLineLength, fp[1]) != NULL) { sumtTokenP = &s[0]; do { GetSumtToken (&tokenType); if (IsSame("[", sumtToken) == SAME) inSumtComment = YES; if (IsSame("]", sumtToken) == SAME) inSumtComment = NO; if (inSumtComment == NO) { if (foundBegin == YES) { if (IsSame("Trees", sumtToken) == SAME) { numTreesInBlock = 0; inTreeBlock = YES; foundBegin = NO; lastTreeBlockBegin = lineNum; } } else { if (IsSame("Begin", sumtToken) == SAME) { if (foundBegin == YES) { MrBayesPrint ("%s Found inappropriate \"Begin\" statement in file\n", spacer); blockErrors = YES; } foundBegin = YES; } else if (IsSame("End", sumtToken) == SAME) { if (inTreeBlock == YES) { numTreeBlocks++; inTreeBlock = NO; lastTreeBlockEnd = lineNum; } else { MrBayesPrint ("%s Found inappropriate \"End\" statement in file\n", spacer); blockErrors = YES; } numTreesInLastBlock = numTreesInBlock; } else if (IsSame("Tree", sumtToken) == SAME) { if (inTreeBlock == YES) numTreesInBlock++; else { MrBayesPrint ("%s Found a \"Tree\" statement that is not in a tree block\n", spacer); blockErrors = YES; } } } } } while (*sumtToken); lineNum++; } /* Now, check some aspects of the tree file, such as the number of tree blocks and whether they are properly terminated. */ if (inTreeBlock == YES) { MrBayesPrint ("%s Unterminated tree block in file %s. You probably need to\n", spacer, comptreeParams.comptFileName2); MrBayesPrint ("%s add a new line to the end of the file with \"End;\" on it.\n", spacer); goto errorExit; } if (inSumtComment == YES) { MrBayesPrint ("%s Unterminated comment in file %s\n", spacer, comptreeParams.comptFileName2); goto errorExit; } if (blockErrors == YES) { MrBayesPrint ("%s Found formatting errors in file %s\n", spacer, comptreeParams.comptFileName2); goto errorExit; } if (lastTreeBlockEnd < lastTreeBlockBegin) { MrBayesPrint ("%s Problem reading tree file %s\n", spacer, comptreeParams.comptFileName2); goto errorExit; } if (numTreesInLastBlock <= 0) { MrBayesPrint ("%s No trees were found in last tree block of file %s\n", spacer, comptreeParams.comptFileName2); goto errorExit; } if (sumtParams.sumtBurnIn > numTreesInLastBlock) { MrBayesPrint ("%s No trees are sampled as the burnin exceeds the number of trees in last block\n", spacer); MrBayesPrint ("%s Try setting burnin to a number less than %d\n", spacer, numTreesInLastBlock); goto errorExit; } /* tell the user that everything is fine */ if (numTreeBlocks == 1) MrBayesPrint ("%s Found one tree block in file \"%s\" with %d trees in last block\n", spacer, comptreeParams.comptFileName2, numTreesInLastBlock); else { MrBayesPrint ("%s Found %d tree blocks in file \"%s\" with %d trees in last block\n", spacer, numTreeBlocks, comptreeParams.comptFileName2, numTreesInLastBlock); MrBayesPrint ("%s Only the %d trees in last tree block will be summarized\n", spacer, numTreesInLastBlock); } /* Now we read the file for real. First, rewind file pointer to beginning of file... */ (void)fseek(fp[1], 0L, 0); /* ...and fast forward to beginning of last tree block. */ for (i=0; i 0) { if (numExcludedTaxa == 1) MrBayesPrint ("%s The following species was absent from trees:\n", spacer); else MrBayesPrint ("%s The following %d species were absent from trees:\n", spacer, numExcludedTaxa); MrBayesPrint ("%s ", spacer); j = lineWidth = 0; for (i=0; i 60) { MrBayesPrint ("\n%s ", spacer); lineWidth = 0; } if (numExcludedTaxa == 1) MrBayesPrint ("%s\n", tempStr); else if (numExcludedTaxa == 2 && j == 1) MrBayesPrint ("%s ", tempStr); else if (j == numExcludedTaxa) MrBayesPrint ("and %s\n", tempStr); else MrBayesPrint ("%s, ", tempStr); } } MrBayesPrint ("\n"); } /* print out information on pruned taxa */ numExcludedTaxa = 0; for (i=0; i 0) { if (numExcludedTaxa == 1) MrBayesPrint ("%s The following species was pruned from trees:\n", spacer); else MrBayesPrint ("%s The following %d species were pruned from trees:\n", spacer, numExcludedTaxa); MrBayesPrint ("%s ", spacer); j = lineWidth = 0; for (i=0; i 60) { MrBayesPrint ("\n%s ", spacer); lineWidth = 0; } if (numExcludedTaxa == 1) MrBayesPrint ("%s\n", tempStr); else if (numExcludedTaxa == 2 && j == 1) MrBayesPrint ("%s ", tempStr); else if (j == numExcludedTaxa) MrBayesPrint ("and %s\n", tempStr); else MrBayesPrint ("%s, ", tempStr); } } MrBayesPrint ("\n"); } /* tell user how many trees were successfully read */ MrBayesPrint ("%s Read %d trees from last tree block (sampling %d of them)\n", spacer, numCompTrees[1], numCompTreesSampled[1]); /* summarize information ***************************************************************************/ if (numIncludedTaxa1 != numIncludedTaxa2) { MrBayesPrint ("%s The number of taxa to be compared in each file is not the same\n", spacer); goto errorExit; } /* how many taxon bipartitions are there for each tree */ if (isSumtTreeRooted == YES) nCompPartitions = 2 * numIncludedTaxa1 - 2; else nCompPartitions = 2 * numIncludedTaxa1 - 3; /* Check that at least one tree was read in. */ if (numSumTreesSampled <= 0) { MrBayesPrint ("%s No trees read in\n", spacer); goto errorExit; } /* Sort partitions... */ if (memAllocs[ALLOC_PARTORIGORDER] == YES) { MrBayesPrint ("%s numFoundOfThisPart not free in AllocBits\n", spacer); goto errorExit; } partOrigOrder = (int *)SafeMalloc((size_t) (numTreePartsFound * sizeof(int))); if (!partOrigOrder) { MrBayesPrint ("%s Problem allocating partOrigOrder (%d)\n", spacer, numTreePartsFound * sizeof(int)); goto errorExit; } memAllocs[ALLOC_PARTORIGORDER] = YES; for (i=0; i= sumtParams.freqDisplay) { MrBayesPrint ("%s %4d -- ", spacer, i+1); ShowParts (&x[0], numIncludedTaxa); MrBayesPrint (" %4d %4d %1.3lf %1.3lf\n", numFoundOfThisPart1[i], numFoundOfThisPart2[i], (MrBFlt)numFoundOfThisPart1[i]/numCompTreesSampled[0], (MrBFlt)numFoundOfThisPart2[i]/numCompTreesSampled[1]); MrBayesPrintf (fpCompParts, "%d\t%d\t%d\t%1.3lf\t%1.3lf\n", i+1, numFoundOfThisPart1[i], numFoundOfThisPart2[i], (MrBFlt)numFoundOfThisPart1[i]/numCompTreesSampled[0], (MrBFlt)numFoundOfThisPart2[i]/numCompTreesSampled[1]); } x += taxonLongsNeeded; } /* make a nifty graph plotting frequencies of clades found in the two tree files */ MrBayesPrint (" \n"); MrBayesPrint ("%s Bivariate plot of clade probabilities: \n", spacer); MrBayesPrint (" \n"); MrBayesPrint ("%s This graph plots the probabilities of clades found in file 1 (the x-axis) \n", spacer); MrBayesPrint ("%s against the probabilities of the same clades found in file 2 (the y-axis). \n", spacer); MrBayesPrint (" \n"); xInc = (MrBFlt) (1.0 / 80.0); yInc = (MrBFlt) (1.0 / 40.0); yUpper = 1.0; yLower = yUpper - yInc; for (yaxis=39; yaxis>=0; yaxis--) { xLower = 0.0; xUpper = xLower + xInc; for (xaxis=0; xaxis<80; xaxis++) { starHolder[xaxis] = 0; for (i=0; i xLower && xProb <= xUpper && yProb > yLower && yProb <= yUpper) starHolder[xaxis] = 1; } xLower += xInc; xUpper = xLower + xInc; } MrBayesPrint ("%s ", spacer); for (xaxis=0; xaxis<80; xaxis++) { prCh = ' '; if ((xaxis == 0 && yaxis == 0) || (xaxis == 79 && yaxis == 39)) prCh = '+'; else if ((xaxis == 0 && yaxis == 39) || (xaxis == 79 && yaxis == 0)) prCh = '+'; else if ((yaxis == 0 || yaxis == 39) && xaxis > 0 && xaxis < 79) prCh = '-'; else if ((xaxis == 0 || xaxis == 79) && yaxis > 0 && yaxis < 39) prCh = '|'; if (starHolder[xaxis] == 1) prCh = '*'; MrBayesPrint ("%c", prCh); } if (yaxis == 39) MrBayesPrint (" 1.00\n"); else if (yaxis == 0) MrBayesPrint (" 0.00\n"); else MrBayesPrint ("\n"); yUpper -= yInc; yLower = yUpper - yInc; } MrBayesPrint (" ^ ^\n"); MrBayesPrint (" 0.00 1.00\n"); /* get tree-to-tree distances */ minNumTrees = numFullCompTreesFound[0]; if (numFullCompTreesFound[1] < minNumTrees) minNumTrees = numFullCompTreesFound[1]; if (memAllocs[ALLOC_TOPO_DISTANCES] == YES) { MrBayesPrint ("%s Topological distances all ready allocated\n", spacer); goto errorExit; } dT1 = (MrBFlt *)SafeMalloc((size_t) (minNumTrees * sizeof(MrBFlt))); if (!dT1) { MrBayesPrint ("%s Problem allocating topological distances\n", spacer); goto errorExit; } dT2 = (MrBFlt *)SafeMalloc((size_t) (minNumTrees * sizeof(MrBFlt))); if (!dT2) { MrBayesPrint ("%s Problem allocating topological distances\n", spacer); goto errorExit; } dT3 = (MrBFlt *)SafeMalloc((size_t) (minNumTrees * sizeof(MrBFlt))); if (!dT3) { MrBayesPrint ("%s Problem allocating topological distances\n", spacer); goto errorExit; } memAllocs[ALLOC_TOPO_DISTANCES] = YES; for (i=0; i maxX) maxX = xVal; if (yVal > maxY) maxY = yVal; } for (i=0; i= screenWidth) k = screenWidth - 1; meanY[k] += yVal; numY[k]++; } MrBayesPrint ("\n +"); for (i=0; i=0; j--) { MrBayesPrint (" |"); for (i=0; i 0) { if (meanY[i] / numY[i] > (((maxY - minY)/screenHeigth)*j)+minY && meanY[i] / numY[i] <= (((maxY - minY)/screenHeigth)*(j+1))+minY) MrBayesPrint ("*"); else MrBayesPrint (" "); } else { MrBayesPrint (" "); } } MrBayesPrint ("|\n"); } MrBayesPrint (" +"); for (i=0; i maxX) maxX = xVal; if (yVal > maxY) maxY = yVal; } for (i=0; i= screenWidth) k = screenWidth - 1; meanY[k] += yVal; numY[k]++; } MrBayesPrint ("\n +"); for (i=0; i=0; j--) { MrBayesPrint (" |"); for (i=0; i 0) { if (meanY[i] / numY[i] > (((maxY - minY)/screenHeigth)*j)+minY && meanY[i] / numY[i] <= (((maxY - minY)/screenHeigth)*(j+1))+minY) MrBayesPrint ("*"); else MrBayesPrint (" "); } else { MrBayesPrint (" "); } } MrBayesPrint ("|\n"); } MrBayesPrint (" +"); for (i=0; i 1) sprintf (fileName,"%s.tree%d", sumtParams.sumtFileName, treeNo+1); else strcpy (fileName, sumtParams.sumtFileName); if (sumtParams.numRuns == 1) MrBayesPrint ("%s Summarizing trees in file \"%s.t\"\n", spacer, fileName); else if (sumtParams.numRuns == 2) MrBayesPrint ("%s Summarizing trees in files \"%s.run1.t\" and \"%s.run2.t\"\n", spacer, fileName, fileName); else if (sumtParams.numRuns > 2) MrBayesPrint ("%s Summarizing trees in files \"%s.run1.t\", \"%s.run2.t\", etc\n", spacer, fileName, fileName); for (runNo=0; runNo < sumtParams.numRuns; runNo++) { /* initialize tree counters */ numSumtTrees = numSumTreesSampled = 0; /* open binary file */ if (sumtParams.numRuns == 1) sprintf (tempName, "%s.t", fileName); else sprintf (tempName, "%s.run%d.t", fileName, runNo+1); if ((fp = OpenBinaryFileR(tempName)) == NULL) { if (strcmp (fileName+strlen(fileName)-2, ".t") == 0) { MrBayesPrint ("%s You probably need to remove '.t' from 'Filename'\n", spacer); MrBayesPrint ("%s Also make sure that 'Nruns' and 'Ntrees' are set correctly\n", spacer); } else MrBayesPrint ("%s Make sure that 'Nruns' and 'Ntrees' are set correctly\n", spacer); goto errorExit; } /* find out what type of line termination is used */ if (runNo == 0 && treeNo == 0) { lineTerm = LineTermType (fp); if (lineTerm == LINETERM_MAC) MrBayesPrint ("%s Macintosh line termination\n", spacer); else if (lineTerm == LINETERM_DOS) MrBayesPrint ("%s DOS line termination\n", spacer); else if (lineTerm == LINETERM_UNIX) MrBayesPrint ("%s UNIX line termination\n", spacer); else { MrBayesPrint ("%s Unknown line termination\n", spacer); goto errorExit; } if (sumtParams.numRuns > 1) MrBayesPrint ("%s Examining first file ...\n", spacer); else MrBayesPrint ("%s Examining file ...\n", spacer); } else if (LineTermType (fp) != lineTerm) { MrBayesPrint ("%s Inconsistent line termination for file %d\n", spacer, runNo + 1); goto errorExit; } /* find length of longest line */ longestLineLength = LongestLine (fp); longestLineLength += 10; /* allocate a string long enough to hold a line */ if (runNo == 0 && treeNo == 0) { if (memAllocs[ALLOC_SUMTSTRING] == YES) { MrBayesPrint ("%s Sumt string is already allocated\n", spacer); goto errorExit; } s = (char *)SafeMalloc((size_t) (longestLineLength * sizeof(char))); if (!s) { MrBayesPrint ("%s Problem allocating string for reading sumt file\n", spacer); goto errorExit; } memAllocs[ALLOC_SUMTSTRING] = YES; } else { free (s); s = (char *) SafeMalloc (sizeof (char) * longestLineLength); } /* close binary file */ SafeFclose (&fp); /* open text file */ if ((fp = OpenTextFileR(tempName)) == NULL) goto errorExit; /* Check file for appropriate blocks. We want to find the last tree block in the file and start from there. */ foundBegin = inTreeBlock = blockErrors = inSumtComment = NO; lineNum = numTreesInBlock = lastTreeBlockBegin = lastTreeBlockEnd = numTreeBlocks = numTreesInLastBlock = 0; while (fgets (s, longestLineLength, fp) != NULL) { sumtTokenP = &s[0]; do { GetSumtToken (&tokenType); if (IsSame("[", sumtToken) == SAME) inSumtComment = YES; if (IsSame("]", sumtToken) == SAME) inSumtComment = NO; if (inSumtComment == NO) { if (foundBegin == YES) { if (IsSame("Trees", sumtToken) == SAME) { numTreesInBlock = 0; inTreeBlock = YES; foundBegin = NO; lastTreeBlockBegin = lineNum; } } else { if (IsSame("Begin", sumtToken) == SAME) { if (foundBegin == YES) { MrBayesPrint ("%s Found inappropriate \"Begin\" statement in file\n", spacer); blockErrors = YES; } foundBegin = YES; } else if (IsSame("End", sumtToken) == SAME) { if (inTreeBlock == YES) { numTreeBlocks++; inTreeBlock = NO; lastTreeBlockEnd = lineNum; } else { MrBayesPrint ("%s Found inappropriate \"End\" statement in file\n", spacer); blockErrors = YES; } numTreesInLastBlock = numTreesInBlock; } else if (IsSame("Tree", sumtToken) == SAME) { if (inTreeBlock == YES) { numTreesInBlock++; if (numTreesInBlock == 1) { sumtBrlensDef = NO; for (i=0; s[i]!='\0'; i++) { if (s[i] == ':') { sumtBrlensDef = YES; break; } } } } else { MrBayesPrint ("%s Found a \"Tree\" statement that is not in a tree block\n", spacer); blockErrors = YES; } } } } } while (*sumtToken); lineNum++; } if (runNo == 0) firstSumtBrlensDef = sumtBrlensDef; else if (firstSumtBrlensDef != sumtBrlensDef) { MrBayesPrint ("%s Tree files with and without brlens mixed\n"); goto errorExit; } /* Now, check some aspects of the tree file, such as the number of tree blocks and whether they are properly terminated. */ if (inTreeBlock == YES) { MrBayesPrint ("%s Unterminated tree block in file %s. You probably need to\n", spacer, tempName); MrBayesPrint ("%s add a new line to the end of the file with \"End;\" on it.\n", spacer); goto errorExit; } if (inSumtComment == YES) { MrBayesPrint ("%s Unterminated comment in file %s\n", spacer, tempName); goto errorExit; } if (blockErrors == YES) { MrBayesPrint ("%s Found formatting errors in file %s\n", spacer, tempName); goto errorExit; } if (lastTreeBlockEnd < lastTreeBlockBegin) { MrBayesPrint ("%s Problem reading tree file %s\n", spacer, tempName); goto errorExit; } if (numTreesInLastBlock <= 0) { MrBayesPrint ("%s No trees were found in last tree block of file %s\n", spacer, tempName); goto errorExit; } if (sumtParams.sumtBurnIn > numTreesInLastBlock) { MrBayesPrint ("%s No trees are sampled as the burnin exceeds the number of trees in last block\n", spacer); MrBayesPrint ("%s Try setting burnin to a number less than %d\n", spacer, numTreesInLastBlock); goto errorExit; } /* tell the user that everything is fine */ if (runNo == 0) { if (numTreeBlocks == 1) MrBayesPrint ("%s Found one tree block in file \"%s\" with %d trees in last block\n", spacer, tempName, numTreesInLastBlock); else { MrBayesPrint ("%s Found %d tree blocks in file \"%s\" with %d trees in last block\n", spacer, numTreeBlocks, tempName, numTreesInLastBlock); MrBayesPrint ("%s Only the %d trees in last tree block will be summarized\n", spacer, numTreesInLastBlock); } if (sumtParams.numRuns > 1) MrBayesPrint ("%s Expecting the same number of trees in the last tree block of all files\n", spacer); firstFileNumTrees = numTreesInLastBlock; } else { if (numTreesInLastBlock != firstFileNumTrees) { MrBayesPrint ("%s Found %d trees in first file but %d trees in file %d\n", spacer, firstFileNumTrees, numTreesInLastBlock); goto errorExit; } } /* Now we read the file for real. First, rewind file pointer to beginning of file... */ (void)fseek(fp, 0L, 0); /* ...and fast forward to beginning of last tree block. */ for (i=0; i 1) MrBayesPrint ("\n%s Tree reading status for tree %d:\n\n", spacer, treeNo+1); else MrBayesPrint ("\n%s Tree reading status:\n\n", spacer); MrBayesPrint ("%s 0 10 20 30 40 50 60 70 80 90 100\n", spacer); MrBayesPrint ("%s v-------v-------v-------v-------v-------v-------v-------v-------v-------v-------v\n", spacer); MrBayesPrint ("%s *", spacer); numAsterices = 0; } /* ...and parse file, tree-by-tree. We are only parsing lines between the "begin trees" and "end" statements. We don't actually get those lines, however, but rather the lines between those statements. */ expecting = Expecting(COMMAND); inSumtBlock = YES; isTranslateDef = NO; numTranslates = whichTranslate = 0; for (i=0; i 1) { for (i=0; i 0) { if (numExcludedTaxa == 1) MrBayesPrint ("%s The following species was absent from trees:\n", spacer); else MrBayesPrint ("%s The following %d species were absent from trees:\n", spacer, numExcludedTaxa); MrBayesPrint ("%s ", spacer); j = lineWidth = 0; for (i=0; i 60) { MrBayesPrint ("\n%s ", spacer); lineWidth = 0; } if (numExcludedTaxa == 1) MrBayesPrint ("%s\n", tempStr); else if (numExcludedTaxa == 2 && j == 1) MrBayesPrint ("%s ", tempStr); else if (j == numExcludedTaxa) MrBayesPrint ("and %s\n", tempStr); else MrBayesPrint ("%s, ", tempStr); } } MrBayesPrint ("\n"); } /* print out information on pruned taxa */ numExcludedTaxa = 0; for (i=0; i 0) { if (numExcludedTaxa == 1) MrBayesPrint ("%s The following species was pruned from trees:\n", spacer); else MrBayesPrint ("%s The following %d species were pruned from trees:\n", spacer, numExcludedTaxa); MrBayesPrint ("%s ", spacer); j = lineWidth = 0; for (i=0; i 60) { MrBayesPrint ("\n%s ", spacer); lineWidth = 0; } if (numExcludedTaxa == 1) MrBayesPrint ("%s\n", tempStr); else if (numExcludedTaxa == 2 && j == 1) MrBayesPrint ("%s ", tempStr); else if (j == numExcludedTaxa) MrBayesPrint ("and %s\n", tempStr); else MrBayesPrint ("%s, ", tempStr); } } MrBayesPrint ("\n"); } } /* Update total number of sampled trees */ numTotalTreesSampled += numSumTreesSampled; /* tell user how many trees were successfully read */ if (sumtParams.numRuns == 1) MrBayesPrint ("%s Read %d trees from last tree block (sampling %d of them)\n", spacer, numSumtTrees, numSumTreesSampled); else if (sumtParams.numRuns > 1) { if (runNo == 0) firstFileNumSumTreesSampled = numSumTreesSampled; else if (numSumTreesSampled != firstFileNumSumTreesSampled) { if (runNo != sumtParams.numRuns - 1) MrBayesPrint ("%s Found %d post-burnin trees in the first file but %d post-burnin trees in file %d\n", spacer, firstFileNumSumTreesSampled, numSumTreesSampled, runNo+1); else MrBayesPrint ("\n\n%s Found %d post-burnin trees in the first file but %d post-burnin trees in file %d\n", spacer, firstFileNumSumTreesSampled, numSumTreesSampled, runNo+1); goto errorExit; } if (runNo == sumtParams.numRuns - 1) { MrBayesPrint ("%s Read a total of %d trees in %d files (sampling %d of them)\n", spacer, sumtParams.numRuns * numSumtTrees, sumtParams.numRuns, numTotalTreesSampled); MrBayesPrint ("%s (Each file contained %d trees of which %d were sampled)\n", spacer, numSumtTrees, numSumTreesSampled); } } /* Check that at least one tree was read in. */ if (numSumTreesSampled <= 0) { MrBayesPrint ("%s No trees read in\n", spacer); goto errorExit; } SafeFclose (&fp); } /* next run for this tree */ /* Sort partitions... */ if (treeNo == 0) { if (memAllocs[ALLOC_PARTORIGORDER] == YES) { MrBayesPrint ("%s partOrigOrder not free in DoSumt\n", spacer); goto errorExit; } partOrigOrder = (int *)SafeMalloc((size_t) (numTreePartsFound * sizeof(int))); if (!partOrigOrder) { MrBayesPrint ("%s Problem allocating partOrigOrder (%d)\n", spacer, numTreePartsFound * sizeof(int)); goto errorExit; } memAllocs[ALLOC_PARTORIGORDER] = YES; } else { free (partOrigOrder); partOrigOrder = (int *)SafeMalloc((size_t) (numTreePartsFound * sizeof(int))); if (!partOrigOrder) { MrBayesPrint ("%s Problem allocating partOrigOrder (%d)\n", spacer, numTreePartsFound * sizeof(int)); goto errorExit; } } for (i=0; i longestName) longestName = len; } MrBayesPrint (" \n"); MrBayesPrint ("%s General explanation: \n", spacer); MrBayesPrint (" \n"); MrBayesPrint ("%s A taxon bibartition is specified by removing a branch, thereby dividing the \n", spacer); MrBayesPrint ("%s species into those to the left and those to the right of the branch. Here, \n", spacer); MrBayesPrint ("%s taxa to one side of the removed branch are denoted \".\" and those to the \n", spacer); MrBayesPrint ("%s other side are denoted \"*\". The output includes the bipartition number \n", spacer); MrBayesPrint ("%s (ID; sorted from highest to lowest probability), bipartition (e.g., ...**..), \n", spacer); MrBayesPrint ("%s number of times the bipartition was observed (#obs), the posterior probabil- \n", spacer); MrBayesPrint ("%s ity of the bipartition, and, if branch lengths were recorded on the trees in \n", spacer); MrBayesPrint ("%s the file, the average (Mean(v)) and variance (Var(v)) of the lengths. Each \n", spacer); MrBayesPrint ("%s \".\" or \"*\" in the bipartition represents a taxon that is to the left or \n", spacer); MrBayesPrint ("%s right of the removed branch. A list of the taxa in the bipartition is given \n", spacer); MrBayesPrint ("%s before the list of bipartitions. If you summarize several independent analy- \n", spacer); MrBayesPrint ("%s ses, convergence diagnostics are presented for both the posterior probabil- \n", spacer); MrBayesPrint ("%s ities of bipartitions (bipartition or split frequencies) and branch lengths \n", spacer); MrBayesPrint ("%s (if recorded on the trees in the files). In the former case, the diagnostic is\n", spacer); MrBayesPrint ("%s the standard deviation of the partition frequencies (Stdev(s)), in the second \n", spacer); MrBayesPrint ("%s case it is the potential scale reduction factor (PSRF) of Gelman and Rubin \n", spacer); MrBayesPrint ("%s (1992). Stdev(s) is expected to approach 0 and PSRF is expected to approach 1 \n", spacer); MrBayesPrint ("%s as runs converge onto the posterior probability distribution. Note that these \n", spacer); MrBayesPrint ("%s values may be unreliable if the partition is not present in all runs (the \n", spacer); MrBayesPrint ("%s last column indicates the number of runs that sampled the partition if more \n", spacer); MrBayesPrint ("%s than one run is summarized). The PSRF is also sensitive to small sample \n", spacer); MrBayesPrint ("%s sizes and it should only be considered a rough guide to convergence since \n", spacer); MrBayesPrint ("%s some of the assumptions allowing one to interpret it as a true potential \n", spacer); MrBayesPrint ("%s scale reduction factor are violated in the phylogenetic context. \n", spacer); MrBayesPrint ("%s \n", spacer); MrBayesPrint ("%s List of taxa in bipartitions: \n", spacer); MrBayesPrint (" \n"); j = 1; for (k=0; k 1) { MrBayesPrint ("%s Results for tree number %d\n", spacer, treeNo+1); MrBayesPrint ("%s ==========================\n\n", spacer); } MrBayesPrint (" \n"); MrBayesPrint ("%s Summary statistics for taxon bipartitions: \n\n", spacer); /* calculate a couple of numbers that are handy to have */ numTreePartsToPrint = 0; for (i=0; i= sumtParams.freqDisplay) numTreePartsToPrint++; } maxWidthID = (int) (log10 (numTreePartsToPrint)) + 1; if (maxWidthID < 2) maxWidthID = 2; maxWidthNumberPartitions = (int) (log10 (numFoundOfThisPart[0])) + 1; if (maxWidthNumberPartitions < 4) maxWidthNumberPartitions = 4; /* print header to screen */ MrBayesPrint ("%s ", spacer); MrBayesPrint ("%*s -- Partition", maxWidthID, "ID"); tableWidth = maxWidthID + 13; for (i=9; i 1) { MrBayesPrint (" Stdev(s)"); tableWidth += 10; } if (sumtBrlensDef == YES) { MrBayesPrint (" Mean(v) Var(v) "); tableWidth += 20; if (sumtParams.numRuns > 1) { MrBayesPrint (" PSRF"); tableWidth += 7; } } if (sumtParams.numRuns > 1) { MrBayesPrint (" Nruns"); tableWidth += 7; } MrBayesPrint ("\n%s ", spacer); for (i=0; i 1) { sum_s = 0.0; sumsq_s = 0.0; for (n=j=0; n 0) j++; f = (MrBFlt) numFoundInRunOfPart[n][i] / (MrBFlt) firstFileNumSumTreesSampled; sum_s += f; sumsq_s += f * f; } var_s = sumsq_s - sum_s * sum_s / (MrBFlt) sumtParams.numRuns; var_s /= (sumtParams.numRuns - 1); if (var_s > 0.0) stddev_s = sqrt (var_s); else stddev_s = 0.0; if (sumtBrlensDef == YES) { if (numFoundOfThisPart[i] - j == 0) varW = 0.0; else varW = sWithinBrlens[i] / (MrBFlt) (numFoundOfThisPart[i] - j); if (j == 1) varB = 0.0; else varB = (sumsqB[i] - sumB[i]*sumB[i]/(MrBFlt) j) / (MrBFlt) (j - 1); if (varW > 0.0) sqrt_R = sqrt ((MrBFlt)(firstFileNumSumTreesSampled - 1) / (MrBFlt) (firstFileNumSumTreesSampled) + ((MrBFlt) (j + 1) / (MrBFlt) (j)) * (varB / varW)); else sqrt_R = -1.0; } if (j == sumtParams.numRuns) unreliable = NO; else { unreliable = YES; oneUnreliable = YES; } } f = (MrBFlt) numFoundOfThisPart[i] / (MrBFlt) numTotalTreesSampled; MrBayesPrint (" %*d %1.6lf", maxWidthNumberPartitions, numFoundOfThisPart[i], f); if (sumtParams.numRuns > 1) MrBayesPrint (" %1.6lf", stddev_s); if (sumtBrlensDef == YES) { MrBayesPrint (" %1.6lf %1.6lf", aBrlens[i], var); if (sumtParams.numRuns > 1) { if (varW <= 0.0) MrBayesPrint (" > 2.0"); else MrBayesPrint (" %1.3lf", sqrt_R); } } if (sumtParams.numRuns > 1) { MrBayesPrint (" %3d", j); } if (unreliable == YES) MrBayesPrint (" *\n"); else MrBayesPrint ("\n"); x += taxonLongsNeeded; } MrBayesPrint ("%s ", spacer); for (i=0; i 1) { MrBayesPrintf (fpParts, "[ STDDEV(s) = Standard deviation of partition probabilities across partitions]\n"); } if (sumtBrlensDef == YES) { MrBayesPrintf (fpParts, "[ BRLEN = Mean branch length]\n"); MrBayesPrintf (fpParts, "[ VAR = Variance of branch length]\n"); if (sumtParams.numRuns > 1) { MrBayesPrintf (fpParts, "[ PSRF = Potential scale reduction factor for branch lengths]\n"); } } if (sumtParams.numRuns > 1) { MrBayesPrintf (fpParts, "[ NRUNS = Number of runs with this partition]\n"); } MrBayesPrintf (fpParts, "\nID\tPART\tNUM\tPROB"); if (sumtParams.numRuns > 1) MrBayesPrintf (fpParts, "\tSTDDEV(s)"); if (sumtBrlensDef == YES) { MrBayesPrintf (fpParts, "\tBRLEN\tVAR"); if (sumtParams.numRuns > 1) MrBayesPrintf (fpParts, "\tR"); } if (sumtParams.numRuns > 1) MrBayesPrintf (fpParts, "NRUNS\n"); else MrBayesPrintf (fpParts, "\n"); /* then print partitions in tab-delimited format */ x = &treePartsFound[0]; for (i=0; i 1) { sum_s = 0.0; sumsq_s = 0.0; for (n=j=0; n 0) j++; f = (MrBFlt) numFoundInRunOfPart[n][i] / (MrBFlt) firstFileNumSumTreesSampled; sum_s += f; sumsq_s += f; } if (j == 1) var_s = 0.0; else var_s = sumsq_s - sum_s * sum_s / (MrBFlt) sumtParams.numRuns; if (j == 1) var_s = 0.0; else var_s = var_s / (sumtParams.numRuns - j); if (var_s > 0.0) stddev_s = sqrt (var_s); else stddev_s = 0.0; if (sumtBrlensDef == YES) { if (numFoundOfThisPart[i] - j == 0) varW = 0.0; else varW = sWithinBrlens[i] / (MrBFlt) (numFoundOfThisPart[i] - j); if (j == 1) varB = 0.0; else varB = (sumsqB[i] - sumB[i]*sumB[i]/(MrBFlt) j) / (MrBFlt) (j - 1); if (varW > 0.0) sqrt_R = sqrt ((MrBFlt)(firstFileNumSumTreesSampled - 1) / (MrBFlt) (firstFileNumSumTreesSampled) + ((MrBFlt) (j + 1) / (MrBFlt) (j)) * (varB / varW)); else sqrt_R = -1.0; } if (j == sumtParams.numRuns) unreliable = NO; else { unreliable = YES; oneUnreliable = YES; } } f = (MrBFlt) numFoundOfThisPart[i] / (MrBFlt) numTotalTreesSampled; MrBayesPrintf (fpParts, "%d\t%1.6lf", numFoundOfThisPart[i], f); if (sumtParams.numRuns > 1) MrBayesPrintf (fpParts, "\t%1.6lf", stddev_s); if (sumtBrlensDef == YES) { MrBayesPrintf (fpParts, "\t%1.6lf\t%1.6lf", aBrlens[i], var); if (sumtParams.numRuns > 1) { if (varW <= 0.0) MrBayesPrintf (fpParts, "\tN/A"); else MrBayesPrintf (fpParts, "\t%1.3lf", sqrt_R); } } if (sumtParams.numRuns > 1) MrBayesPrintf (fpParts, "\t%d\n", j); else MrBayesPrintf (fpParts, "\n"); x += taxonLongsNeeded; } /* get branch lengths and print to file if appropriate */ if (sumtBrlensDef == YES && sumtParams.printBrlensToFile == YES) { if (BrlensVals (treeNo, s, longestLineLength, lastTreeBlockBegin, lastTreeBlockEnd) == ERROR) goto errorExit; } /* make the majority rule consensus tree */ if (ConTree () == ERROR) goto errorExit; /* get probabilities of individual trees */ if (TreeProb () == ERROR) goto errorExit; SafeFclose (&fpParts); SafeFclose (&fpCon); SafeFclose (&fpTrees); SafeFclose (&fpBrlens); } /* next tree */ /* free memory and file pointers */ if (memAllocs[ALLOC_SUMTSTRING] == YES) { free (s); memAllocs[ALLOC_SUMTSTRING] = NO; } if (memAllocs[ALLOC_PARTORIGORDER] == YES) { free (partOrigOrder); memAllocs[ALLOC_PARTORIGORDER] = NO; } FreeBits (); expecting = Expecting(COMMAND); # if defined (MPI_ENABLED) } # endif return (NO_ERROR); /* error exit */ errorExit: expecting = Expecting(COMMAND); if (memAllocs[ALLOC_SUMTSTRING] == YES) { free (s); memAllocs[ALLOC_SUMTSTRING] = NO; } if (memAllocs[ALLOC_PARTORIGORDER] == YES) { free (partOrigOrder); memAllocs[ALLOC_PARTORIGORDER] = NO; } FreeBits (); SafeFclose (&fp); SafeFclose (&fpParts); SafeFclose (&fpCon); SafeFclose (&fpTrees); SafeFclose (&fpBrlens); strcpy (spacer, ""); strcpy (sumtToken, "Sumt"); i = 0; if (FindValidCommand (sumtToken, &i) == ERROR) MrBayesPrint ("%s Could not find sumt\n", spacer); return (ERROR); } int DoSumtParm (char *parmName, char *tkn) { int tempI; MrBFlt tempD; char tempStr[100]; if (defMatrix == NO) { MrBayesPrint ("%s A matrix must be specified before sumt can be used\n", spacer); return (ERROR); } if (expecting == Expecting(PARAMETER)) { expecting = Expecting(EQUALSIGN); } else { if (!strcmp(parmName, "Xxxxxxxxxx")) { expecting = Expecting(PARAMETER); expecting |= Expecting(SEMICOLON); } /* set Filename (sumtParams.sumtFileName) ***************************************************/ else if (!strcmp(parmName, "Filename")) { if (expecting == Expecting(EQUALSIGN)) { expecting = Expecting(ALPHA); readWord = YES; } else if (expecting == Expecting(ALPHA)) { strcpy (sumtParams.sumtFileName, tkn); MrBayesPrint ("%s Setting sumt filename to %s\n", spacer, sumtParams.sumtFileName); expecting = Expecting(PARAMETER) | Expecting(SEMICOLON); } else return (ERROR); } /* set Displaygeq (sumtParams.freqDisplay) *******************************************************/ else if (!strcmp(parmName, "Displaygeq")) { if (expecting == Expecting(EQUALSIGN)) expecting = Expecting(NUMBER); else if (expecting == Expecting(NUMBER)) { sscanf (tkn, "%lf", &tempD); sumtParams.freqDisplay = tempD; MrBayesPrint ("%s Showing partitions with probability greater than or equal to %lf\n", spacer, sumtParams.freqDisplay); expecting = Expecting(PARAMETER) | Expecting(SEMICOLON); } else return (ERROR); } /* set Burnin (sumtParams.sumtBurnIn) *******************************************************/ else if (!strcmp(parmName, "Burnin")) { if (expecting == Expecting(EQUALSIGN)) expecting = Expecting(NUMBER); else if (expecting == Expecting(NUMBER)) { sscanf (tkn, "%d", &tempI); sumtParams.sumtBurnIn = tempI; MrBayesPrint ("%s Setting sumt burnin to %ld\n", spacer, sumtParams.sumtBurnIn); expecting = Expecting(PARAMETER) | Expecting(SEMICOLON); } else return (ERROR); } /* set Nruns (sumtParams.numRuns) *******************************************************/ else if (!strcmp(parmName, "Nruns")) { if (expecting == Expecting(EQUALSIGN)) expecting = Expecting(NUMBER); else if (expecting == Expecting(NUMBER)) { sscanf (tkn, "%d", &tempI); if (tempI < 1) { MrBayesPrint ("%s Nruns must be at least 1\n", spacer); return (ERROR); } else { sumtParams.numRuns = tempI; MrBayesPrint ("%s Setting sumt nruns to %ld\n", spacer, sumtParams.numRuns); expecting = Expecting(PARAMETER) | Expecting(SEMICOLON); } } else return (ERROR); } /* set Ntrees (sumtParams.numTrees) *******************************************************/ else if (!strcmp(parmName, "Ntrees")) { if (expecting == Expecting(EQUALSIGN)) expecting = Expecting(NUMBER); else if (expecting == Expecting(NUMBER)) { sscanf (tkn, "%d", &tempI); if (tempI < 1) { MrBayesPrint ("%s Ntrees must be at least 1\n", spacer); return (ERROR); } else { sumtParams.numTrees = tempI; MrBayesPrint ("%s Setting sumt ntrees to %ld\n", spacer, sumtParams.numTrees); expecting = Expecting(PARAMETER) | Expecting(SEMICOLON); } } else return (ERROR); } /* set Contype (sumtParams.sumtConType) *****************************************************/ else if (!strcmp(parmName, "Contype")) { if (expecting == Expecting(EQUALSIGN)) expecting = Expecting(ALPHA); else if (expecting == Expecting(ALPHA)) { if (IsArgValid(tkn, tempStr) == NO_ERROR) { strcpy (sumtParams.sumtConType, tempStr); MrBayesPrint ("%s Setting sumt contype to %s\n", spacer, sumtParams.sumtConType); } expecting = Expecting(PARAMETER) | Expecting(SEMICOLON); } else return (ERROR); } /* set Calctrprobs (sumtParams.calcTrprobs) *********************************************/ else if (!strcmp(parmName, "Calctrprobs")) { if (expecting == Expecting(EQUALSIGN)) expecting = Expecting(ALPHA); else if (expecting == Expecting(ALPHA)) { if (IsArgValid(tkn, tempStr) == NO_ERROR) { if (!strcmp(tempStr, "Yes")) sumtParams.calcTrprobs = YES; else sumtParams.calcTrprobs = NO; } else { MrBayesPrint ("%s Invalid argument for calctrprobs\n", spacer); return (ERROR); } if (sumtParams.calcTrprobs == YES) MrBayesPrint ("%s Setting calctrprobs to yes\n", spacer); else MrBayesPrint ("%s Setting calctrprobs to no\n", spacer); expecting = Expecting(PARAMETER) | Expecting(SEMICOLON); } else return (ERROR); } /* set Showtreeprobs (sumtParams.showSumtTrees) *********************************************/ else if (!strcmp(parmName, "Showtreeprobs")) { if (expecting == Expecting(EQUALSIGN)) expecting = Expecting(ALPHA); else if (expecting == Expecting(ALPHA)) { if (IsArgValid(tkn, tempStr) == NO_ERROR) { if (!strcmp(tempStr, "Yes")) sumtParams.showSumtTrees = YES; else sumtParams.showSumtTrees = NO; } else { MrBayesPrint ("%s Invalid argument for showtreeprobs\n", spacer); return (ERROR); } if (sumtParams.showSumtTrees == YES) MrBayesPrint ("%s Setting showtreeprobs to yes\n", spacer); else MrBayesPrint ("%s Setting showtreeprobs to no\n", spacer); expecting = Expecting(PARAMETER) | Expecting(SEMICOLON); } else return (ERROR); } /* set Printbrlens (sumtParams.printBrlensToFile) *********************************************/ else if (!strcmp(parmName, "Printbrlens")) { if (expecting == Expecting(EQUALSIGN)) expecting = Expecting(ALPHA); else if (expecting == Expecting(ALPHA)) { if (IsArgValid(tkn, tempStr) == NO_ERROR) { if (!strcmp(tempStr, "Yes")) sumtParams.printBrlensToFile = YES; else sumtParams.printBrlensToFile = NO; } else { MrBayesPrint ("%s Invalid argument for printbrlens\n", spacer); return (ERROR); } if (sumtParams.printBrlensToFile == YES) MrBayesPrint ("%s Setting printbrlens to yes\n", spacer); else MrBayesPrint ("%s Setting printbrlens to no\n", spacer); expecting = Expecting(PARAMETER) | Expecting(SEMICOLON); } else return (ERROR); } /* set Brlensgeq (sumtParams.brlensFreqDisplay) *******************************************************/ else if (!strcmp(parmName, "Brlensgeq")) { if (expecting == Expecting(EQUALSIGN)) expecting = Expecting(NUMBER); else if (expecting == Expecting(NUMBER)) { sscanf (tkn, "%lf", &tempD); sumtParams.brlensFreqDisplay = tempD; MrBayesPrint ("%s Printing branch lengths to file for partitions with probability >= %lf\n", spacer, sumtParams.brlensFreqDisplay); expecting = Expecting(PARAMETER) | Expecting(SEMICOLON); } else return (ERROR); } /* set Ordertaxa (sumtParams.orderTaxa) *********************************************/ else if (!strcmp(parmName, "Ordertaxa")) { if (expecting == Expecting(EQUALSIGN)) expecting = Expecting(ALPHA); else if (expecting == Expecting(ALPHA)) { if (IsArgValid(tkn, tempStr) == NO_ERROR) { if (!strcmp(tempStr, "Yes")) sumtParams.orderTaxa = YES; else sumtParams.orderTaxa = NO; } else { MrBayesPrint ("%s Invalid argument for ordertaxa\n", spacer); return (ERROR); } if (sumtParams.orderTaxa == YES) MrBayesPrint ("%s Setting ordertaxa to yes\n", spacer); else MrBayesPrint ("%s Setting ordertaxa to no\n", spacer); expecting = Expecting(PARAMETER) | Expecting(SEMICOLON); } else return (ERROR); } else return (ERROR); } return (NO_ERROR); } int DoTranslate (void) { if (inSumtBlock == NO) { MrBayesPrint ("%s You must be in a trees block to read a translate command\n", spacer); return (ERROR); } numTranslates++; isTranslateDef = YES; # if 0 MrBayesPrint ("%s Defining a translation table\n", spacer); /* print out translate table */ { int i, j, len, longestLen; char tempName[100]; /*MrBayesPrint ("%s\n", transFrom); MrBayesPrint ("%s\n", transTo);*/ longestLen = 0; for (i=0; i longestLen) longestLen = len; } for (i=0; i numTaxa) { MrBayesPrint ("%s Too many taxa in tree\n", spacer); goto errorExit; } /* Increment number of trees read in. */ numSumtTrees++; if (comparingFiles == YES) numCompTrees[fileNum]++; /* update status bar */ if (comparingFiles == YES) numTreesInThisFile = numCompTrees[fileNum]; else numTreesInThisFile = numSumtTrees; if (numTreesInLastBlock * sumtParams.numRuns < 80) { printEvery = 1; nAstPerPrint = 80 / (numTreesInLastBlock * sumtParams.numRuns); if (numSumtTrees % printEvery == 0) { for (i=0; i sumtParams.sumtBurnIn) || (comparingFiles == YES && numCompTrees[fileNum] > comptreeParams.comptBurnIn)) { /* Increment the number of trees sampled. */ numSumTreesSampled++; if (comparingFiles == YES) numCompTreesSampled[fileNum]++; /* Add a node to the root of the tree if this is a rooted tree. */ if (isSumtTreeRooted == YES) { if (pSumtPtr->anc == NULL) { qSumtPtr = &sumtNodes[nextAvailableSumtNode]; nextAvailableSumtNode++; qSumtPtr->left = pSumtPtr; pSumtPtr->anc = qSumtPtr; pSumtPtr = qSumtPtr; sumtRoot = pSumtPtr; } else { MrBayesPrint ("%s Tree is not rooted correctly\n", spacer); goto errorExit; } } /* Check to see if all of the species in the character matrix (already read in) are present in the trees. We only do this for the first tree read in. */ if (numSumTreesSampled == 1) { if (CheckSumtSpecies () == ERROR) goto errorExit; } /* Prune deleted taxa from list of trees, if necessary */ if (PruneSumt () == ERROR) goto errorExit; if (numSumTreesSampled == 1) { for (i=0; i= 2*numTaxa) { MrBayesPrint ("%s Too many nodes on sumt tree\n", spacer); goto errorExit; } if (pSumtPtr->left == NULL) { pSumtPtr = &sumtNodes[nextAvailableSumtNode]; nextAvailableSumtNode++; qSumtPtr->left = pSumtPtr; pSumtPtr->anc = qSumtPtr; qSumtPtr = pSumtPtr; } else if (pSumtPtr->right == NULL) { pSumtPtr = &sumtNodes[nextAvailableSumtNode]; nextAvailableSumtNode++; qSumtPtr->right = pSumtPtr; pSumtPtr->anc = qSumtPtr; qSumtPtr = pSumtPtr; } else if (pSumtPtr->anc == NULL) { pSumtPtr = &sumtNodes[nextAvailableSumtNode]; nextAvailableSumtNode++; qSumtPtr->anc = pSumtPtr; pSumtPtr->left = qSumtPtr; qSumtPtr = pSumtPtr; sumtRoot = pSumtPtr; pSumtPtr->marked = YES; isSumtTreeRooted = NO; } else { MrBayesPrint ("\n ERROR: Tree is not bifurcating\n"); goto errorExit; } } expecting = Expecting(ALPHA); expecting |= Expecting(NUMBER); expecting |= Expecting(LEFTPAR); } else if (expecting == Expecting(ALPHA)) { if (nextAvailableSumtNode+1 >= 2*numTaxa) { MrBayesPrint ("%s Too many nodes on sumt tree\n", spacer); return (ERROR); } if (isTranslateDef == YES) { /* we are using the translation table */ if (CheckString (tkn, transTo, &howMany) == ERROR) { MrBayesPrint ("%s Could not find taxon %s in list of translation taxa\n", spacer, tkn); goto errorExit; } else { if (GetNameFromString (transFrom, tempName, howMany) == ERROR) { MrBayesPrint ("%s Error getting taxon names \n", spacer); goto errorExit; } else { if (CheckString (tempName, taxaNames, &howMany) == ERROR) { MrBayesPrint ("%s Could not find taxon %s in list of taxa\n", spacer, tkn); goto errorExit; } if (tempSet[howMany-1] == YES) { MrBayesPrint ("%s Taxon name %s already used in tree\n", spacer, tkn); goto errorExit; } else tempSet[howMany-1] = YES; } } } else { /* Check to see if the name is in the list of taxon names. */ if (CheckString (tkn, taxaNames, &howMany) == ERROR) { MrBayesPrint ("%s Could not find taxon %s in list of taxa\n", spacer, tkn); goto errorExit; } if (tempSet[howMany-1] == YES) { MrBayesPrint ("%s Taxon name %s already used in tree\n", spacer, tkn); goto errorExit; } else tempSet[howMany-1] = YES; } if (pSumtPtr->left == NULL) { pSumtPtr = &sumtNodes[nextAvailableSumtNode]; strcpy (pSumtPtr->label, tkn); pSumtPtr->taxonName = YES; pSumtPtr->index = howMany - 1; nextAvailableSumtNode++; qSumtPtr->left = pSumtPtr; pSumtPtr->anc = qSumtPtr; qSumtPtr = pSumtPtr; } else if (pSumtPtr->right == NULL) { pSumtPtr = &sumtNodes[nextAvailableSumtNode]; strcpy (pSumtPtr->label, tkn); pSumtPtr->taxonName = YES; pSumtPtr->index = howMany - 1; nextAvailableSumtNode++; qSumtPtr->right = pSumtPtr; pSumtPtr->anc = qSumtPtr; qSumtPtr = pSumtPtr; } else if (pSumtPtr->anc == NULL) { pSumtPtr = &sumtNodes[nextAvailableSumtNode]; strcpy (pSumtPtr->label, tkn); pSumtPtr->taxonName = YES; pSumtPtr->index = howMany - 1; nextAvailableSumtNode++; qSumtPtr->anc = pSumtPtr; pSumtPtr->left = qSumtPtr; qSumtPtr = pSumtPtr; sumtRoot = pSumtPtr; pSumtPtr->marked = YES; isSumtTreeRooted = NO; } else { MrBayesPrint ("%s Tree is not bifurcating\n", spacer); goto errorExit; } numSumtTaxa++; expecting = Expecting(COMMA); if (sumtBrlensDef == YES) expecting |= Expecting(COLON); expecting |= Expecting(RIGHTPAR); } else if (expecting == Expecting(RIGHTPAR)) { if (pSumtPtr->marked == NO) { if (pSumtPtr->anc != NULL) { pSumtPtr = pSumtPtr->anc; qSumtPtr = pSumtPtr; } else { MrBayesPrint ("%s Cannot go down\n", spacer); goto errorExit; } } else { if (pSumtPtr->left != NULL) { pSumtPtr = pSumtPtr->left; qSumtPtr = pSumtPtr; } else { MrBayesPrint ("%s Cannot go down\n", spacer); goto errorExit; } } expecting = Expecting(COMMA); if (sumtBrlensDef == YES) expecting |= Expecting(COLON); expecting |= Expecting(RIGHTPAR); expecting |= Expecting(SEMICOLON); } else if (expecting == Expecting(COLON)) { foundSumtColon = YES; expecting = Expecting(NUMBER); } else if (expecting == Expecting(COMMA)) { if (pSumtPtr->marked == NO) { if (pSumtPtr->anc != NULL) { pSumtPtr = pSumtPtr->anc; qSumtPtr = pSumtPtr; } else { MrBayesPrint ("%s Cannot go down\n", spacer); goto errorExit; } } else { if (pSumtPtr->left != NULL) { pSumtPtr = pSumtPtr->left; qSumtPtr = pSumtPtr; } else { MrBayesPrint ("%s Cannot go down\n", spacer); goto errorExit; } } expecting = Expecting(ALPHA); expecting |= Expecting(NUMBER); expecting |= Expecting(LEFTPAR); } else if (expecting == Expecting(NUMBER)) { if (foundSumtColon == YES) { /* branch length */ sscanf (tkn, "%lf", &tempD); if (pSumtPtr->marked == NO) pSumtPtr->length = tempD; else { if (pSumtPtr->left != NULL) pSumtPtr->left->length = tempD; else { MrBayesPrint ("%s Cannot assign branch length to left node\n", spacer); goto errorExit; } } foundSumtColon = NO; expecting = Expecting(COMMA); expecting |= Expecting(RIGHTPAR); } else { if (isTranslateDef == YES) { /* we are using the translation table */ if (CheckString (tkn, transTo, &howMany) == ERROR) { MrBayesPrint ("%s Could not find taxon %s in list of translation taxa\n", spacer, tkn); goto errorExit; } else { if (GetNameFromString (transFrom, tempName, howMany) == ERROR) { MrBayesPrint ("%s Error getting partition names \n", spacer); goto errorExit; } else { if (CheckString (tempName, taxaNames, &howMany) == ERROR) { MrBayesPrint ("%s Could not find taxon %s in list of taxa\n", spacer, tkn); goto errorExit; } if (tempSet[howMany-1] == YES) { MrBayesPrint ("%s Taxon name %s already used in tree\n", spacer, tkn); goto errorExit; } else tempSet[howMany-1] = YES; } } tempInt = howMany - 1; } else { /* simply use taxon number */ sscanf (tkn, "%d", &tempInt); if (nextAvailableSumtNode+1 >= 2 * numTaxa) { MrBayesPrint ("%s Too many nodes on sumt tree\n", spacer); goto errorExit; } /* Check to see if the name is in the list of taxon names. */ if (CheckString (tkn, taxaNames, &howMany) == ERROR) { /* The number could not be found as a taxon name in the list of taxon names. We will assume that the user has then input taxa as numbers and not the names. */ if (tempSet[tempInt-1] == YES) { MrBayesPrint ("%s Taxon name %d has already been used in tree\n", spacer, tempInt); goto errorExit; } else tempSet[tempInt-1] = YES; tempInt--; } else { /* The taxon name is in the list of taxon names */ howMany--; if (howMany < 0 || howMany >= numTaxa) { MrBayesPrint ("%s Taxon number is out of range\n", spacer); goto errorExit; } if (tempSet[howMany] == YES) { MrBayesPrint ("%s Taxon %d has already been used in tree\n", spacer, howMany+1); goto errorExit; } else tempSet[howMany] = YES; tempInt = howMany; } if (GetNameFromString (taxaNames, tempName, tempInt+1) == ERROR) { MrBayesPrint ("%s Error getting partition names \n", spacer); goto errorExit; } } if (pSumtPtr->left == NULL) { pSumtPtr = &sumtNodes[nextAvailableSumtNode]; strcpy (pSumtPtr->label, tempName); pSumtPtr->taxonName = YES; pSumtPtr->index = tempInt; nextAvailableSumtNode++; qSumtPtr->left = pSumtPtr; pSumtPtr->anc = qSumtPtr; qSumtPtr = pSumtPtr; } else if (pSumtPtr->right == NULL) { pSumtPtr = &sumtNodes[nextAvailableSumtNode]; strcpy (pSumtPtr->label, tempName); pSumtPtr->taxonName = YES; pSumtPtr->index = tempInt; nextAvailableSumtNode++; qSumtPtr->right = pSumtPtr; pSumtPtr->anc = qSumtPtr; qSumtPtr = pSumtPtr; } else if (pSumtPtr->anc == NULL) { pSumtPtr = &sumtNodes[nextAvailableSumtNode]; strcpy (pSumtPtr->label, tempName); pSumtPtr->taxonName = YES; pSumtPtr->index = tempInt; nextAvailableSumtNode++; qSumtPtr->anc = pSumtPtr; pSumtPtr->left = qSumtPtr; qSumtPtr = pSumtPtr; sumtRoot = pSumtPtr; pSumtPtr->marked = YES; isSumtTreeRooted = NO; } else { MrBayesPrint ("%s Tree is not bifurcating\n", spacer); return (ERROR); } numSumtTaxa++; expecting = Expecting(COMMA); expecting |= Expecting(COLON); expecting |= Expecting(RIGHTPAR); } } return (NO_ERROR); errorExit: inSumtBlock = NO; return (ERROR); MrBayesPrint ("%s", parmName); /* just because I am tired of seeing the unused parameter error msg */ } int FindParts (void) { int i, j, nNodes, partNum; MrBFlt bl; SumtNode **downPass, *p; /* allocate memory for downpass */ downPass = (SumtNode **)SafeMalloc((size_t) (2 * numTaxa * sizeof(SumtNode *))); if (!downPass) { MrBayesPrint ("%s Could not allocate downPass\n", spacer); goto errorExit; } /* get the downpass sequence */ i = 0; GetSumtDownPass (sumtRoot, downPass, &i); nNodes = i; j = 0; for (i=0; ianc != NULL) { if (sumtBrlensDef == YES) bl = p->length; else bl = -1.0; if (p->anc->anc != NULL) { if (PartFinder (&treeBits[(p->index) * taxonLongsNeeded], bl, &partNum) == ERROR) goto errorExit; treePartNums[j] = partNum; treePartLengths[j] = bl; j++; } else if (p->anc->anc == NULL && isSumtTreeRooted == NO) { if (PartFinder (&treeBits[(p->index) * taxonLongsNeeded], bl, &partNum) == ERROR) goto errorExit; treePartNums[j] = partNum; treePartLengths[j] = bl; j++; } } } /* count and sort all of the tree partitions IDs */ numTreeParts = j; SortIndParts (&treePartNums[0], &treePartLengths[0], numTreeParts, YES); free (downPass); # if 0 MrBayesPrint ("Partition IDs (%d): \n", numTreeParts); for (i=0; i numFullTreesAllocated) { numFullTreesAllocated += 500; if (RealloateFullTrees () == ERROR) return (ERROR); } x = &fullTreePartIds[numFullTreesFound * 2 * numTaxa]; for (i=0; ileft, i, isThisTreeRooted); FinishSumtTree (p->right, i, isThisTreeRooted); p->marked = NO; if (p->left == NULL && p->right == NULL && p->anc != NULL) { } else if (p->left != NULL && p->right == NULL && p->anc == NULL) { if (isThisTreeRooted == YES) p->index = (*i)++; } else { p->index = (*i)++; } } } int FirstTaxonInPartition (long *partition, int length) { int i, j, nBits, taxon; long x; nBits = sizeof(long) * 8; for (i=taxon=0; ileft != NULL) { for (q=p->left; q!=NULL; q=q->sib) GetConDownPass(downPass, q, i); } downPass[(*i)++] = p; } /* get the actual down pass sequences */ void GetSumtDownPass (SumtNode *p, SumtNode **dp, int *i) { if (p != NULL ) { GetSumtDownPass (p->left, dp, i); GetSumtDownPass (p->right, dp, i); if (p->left != NULL && p->right != NULL && p->anc != NULL) { dp[(*i)++] = p; } else if (p->left == NULL && p->right == NULL && p->anc != NULL) { dp[(*i)++] = p; } else if (p->left != NULL && p->right == NULL && p->anc == NULL) { dp[(*i)++] = p; } } } int GetPartitions (void) { int i, nNodes; SumtNode **downPass, *p; /* allocate memory for downpass */ downPass = (SumtNode **)SafeMalloc((size_t) (2 * numTaxa * sizeof(SumtNode *))); if (!downPass) { MrBayesPrint ("%s Could not allocate downPass\n", spacer); goto errorExit; } /* set all partitions to 0 */ for (i=0; i<2*numTaxa*taxonLongsNeeded; i++) treeBits[i] = 0; /* get the downpass sequence */ i = 0; GetSumtDownPass (sumtRoot, downPass, &i); nNodes = i; /* find which taxa are in the tree */ for (i=0; ileft == NULL && p->right == NULL) { AssignTipPart (p->index, &treeBits[(p->index) * taxonLongsNeeded]); } else if (p->anc != NULL) { if (p->anc->anc != NULL) AssignIntPart (&treeBits[(p->left->index) * taxonLongsNeeded], &treeBits[(p->right->index) * taxonLongsNeeded], &treeBits[(p->index) * taxonLongsNeeded]); else if (p->anc->anc == NULL && isSumtTreeRooted == NO) AssignIntPart (&treeBits[(p->left->index) * taxonLongsNeeded], &treeBits[(p->right->index) * taxonLongsNeeded], &treeBits[(p->index) * taxonLongsNeeded]); } } # if 0 MrBayesPrint ("Partitions:\n"); for (i=0; ileft == NULL && p->right == NULL) { ShowParts (&treeBits[(p->index) * taxonLongsNeeded], numTaxa); MrBayesPrint ("\n"); } else if (p->anc != NULL) { if (p->anc->anc != NULL) { ShowParts (&treeBits[(p->index) * taxonLongsNeeded], numTaxa); MrBayesPrint ("\n"); } else if (p->anc->anc == NULL && isSumtTreeRooted == NO) { ShowParts (&treeBits[(p->index) * taxonLongsNeeded], numTaxa); MrBayesPrint ("\n"); } } } # endif free (downPass); return (NO_ERROR); errorExit: if (downPass) free (downPass); return (ERROR); } void GetSumtToken (int *tokenType) { int allNumbers; register char *temp; (*tokenType) = 0; temp = sumtToken; while (IsWhite(*sumtTokenP) == 1 || IsWhite(*sumtTokenP) == 2) { if (IsWhite(*sumtTokenP) == 2) { *tokenType = RETURNSYMBOL; foundNewLine = YES; /* MrBayesPrint ("RETURN\n"); */ } ++sumtTokenP; } *tokenType = UNKNOWN_TOKEN_TYPE; if (IsIn(*sumtTokenP,"=")) { *temp++ = *sumtTokenP++; *tokenType = EQUALSIGN; } else if (IsIn(*sumtTokenP,";")) { *temp++ = *sumtTokenP++; *tokenType = SEMICOLON; } else if (IsIn(*sumtTokenP,":")) { *temp++ = *sumtTokenP++; *tokenType = COLON; } else if (IsIn(*sumtTokenP,",")) { *temp++ = *sumtTokenP++; *tokenType = COMMA; } else if (IsIn(*sumtTokenP,"#")) { *temp++ = *sumtTokenP++; *tokenType = POUNDSIGN; } else if (IsIn(*sumtTokenP,"(")) { *temp++ = *sumtTokenP++; *tokenType = LEFTPAR; } else if (IsIn(*sumtTokenP,")")) { *temp++ = *sumtTokenP++; *tokenType = RIGHTPAR; } else if (IsIn(*sumtTokenP,"{")) { *temp++ = *sumtTokenP++; *tokenType = LEFTCURL; } else if (IsIn(*sumtTokenP,"}")) { *temp++ = *sumtTokenP++; *tokenType = RIGHTCURL; } else if (IsIn(*sumtTokenP,"[")) { *temp++ = *sumtTokenP++; *tokenType = LEFTCOMMENT; } else if (IsIn(*sumtTokenP,"]")) { *temp++ = *sumtTokenP++; *tokenType = RIGHTCOMMENT; } else if (IsIn(*sumtTokenP,"?")) { *temp++ = *sumtTokenP++; *tokenType = QUESTIONMARK; } else if (IsIn(*sumtTokenP,"-")) { *temp++ = *sumtTokenP++; *tokenType = DASH; } else if (IsIn(*sumtTokenP,"abcdefghijklmnopqrstuvwxyzABCDEFGHIJKLMNOPQRSTUVWXYZ_0123456789.")) { allNumbers = TRUE; if (!IsIn(*sumtTokenP,"0123456789.")) allNumbers = FALSE; *temp++ = *sumtTokenP++; while(IsIn(*sumtTokenP,"abcdefghijklmnopqrstuvwxyzABCDEFGHIJKLMNOPQRSTUVWXYZ_0123456789.")) { if (!IsIn(*sumtTokenP,"0123456789.")) allNumbers = FALSE; *temp++ = *sumtTokenP++; } if (allNumbers == TRUE) *tokenType = NUMBER; else *tokenType = ALPHA; } else if (IsIn(*sumtTokenP,"*")) { *temp++ = *sumtTokenP++; *tokenType = ASTERISK; } else if (IsIn(*sumtTokenP,"/")) { *temp++ = *sumtTokenP++; *tokenType = FORWARDSLASH; } else if (IsIn(*sumtTokenP,"'\\'")) { *temp++ = *sumtTokenP++; *tokenType = BACKSLASH; } else if (IsIn(*sumtTokenP,"!")) { *temp++ = *sumtTokenP++; *tokenType = EXCLAMATIONMARK; } else if (IsIn(*sumtTokenP,"%")) { *temp++ = *sumtTokenP++; *tokenType = PERCENT; } else if (IsIn(*sumtTokenP,"\"")) { *temp++ = *sumtTokenP++; *tokenType = QUOTATIONMARK; } else if (IsIn(*sumtTokenP,"&~+^$@|{}`><")) { *temp++ = *sumtTokenP++; *tokenType = WEIRD; } *temp = '\0'; } int IsPartCompatible (long *smaller, long *larger, int length) { int i; for (i=0; i 1) sprintf (fileName, "%s.tree%d.brlens", sumtParams.sumtFileName, treeNo+1); else sprintf (fileName, "%s.brlens", sumtParams.sumtFileName); /* open file checking for over-write as appropriate */ if ((fpBrlens = OpenNewMBPrintFile(fileName)) == NULL) return ERROR; /* print unique identifier to the file */ len = (int) strlen (stamp); if (len <= 1) { MrBayesPrintf (fpBrlens, "[ID: None Available]\n"); } else { MrBayesPrintf (fpBrlens, "[ID: %s]\n", stamp); } /* print header */ for (i=0; i 1) { sprintf (pFilename, "%s.tree%d.parts", sumtParams.sumtFileName, treeNo+1); sprintf (cFilename, "%s.tree%d.con", sumtParams.sumtFileName, treeNo+1); sprintf (tFilename, "%s.tree%d.trprobs", sumtParams.sumtFileName, treeNo+1); } else { sprintf (pFilename, "%s.parts", sumtParams.sumtFileName); sprintf (cFilename, "%s.con", sumtParams.sumtFileName); sprintf (tFilename, "%s.trprobs", sumtParams.sumtFileName); } /* Open files checking for over-write as appropriate */ if ((fpParts = OpenNewMBPrintFile(pFilename)) == NULL) return ERROR; if ((fpCon = OpenNewMBPrintFile(cFilename)) == NULL) { SafeFclose (&fpParts); return ERROR; } if (sumtParams.calcTrprobs == YES) { if ((fpTrees = OpenNewMBPrintFile(tFilename)) == NULL) { SafeFclose (&fpParts); SafeFclose (&fpCon); return ERROR; } } /* print #NEXUS if appropriate */ fprintf (fpCon, "#NEXUS\n\n"); if (sumtParams.calcTrprobs == YES) fprintf (fpTrees, "#NEXUS\n\n"); /* print unique identifiers to each file */ len = (int) strlen (stamp); if (len <= 1) { if (sumtParams.calcTrprobs == YES) fprintf (fpParts, "[ID: None Available]\n"); fprintf (fpCon, "[ID: None Available]\n"); fprintf (fpTrees, "[ID: None Available]\n"); } else { if (sumtParams.calcTrprobs == YES) fprintf (fpParts, "[ID: %s]\n", stamp); fprintf (fpCon, "[ID: %s]\n", stamp); fprintf (fpTrees, "[ID: %s]\n", stamp); } return (NO_ERROR); } int PartFinder (long *p, MrBFlt bl, int *partID) { int i, n, foundPart, nDiff, whichPart; long *x; MrBFlt f; if (printingBrlens == YES) { foundPart = NO; x = &treePartsFound[0]; for (n=0; n 1) numFoundInRunOfPart[runIndex][n]++; if (comparingFiles == YES) { if (fileNum == 0) numFoundOfThisPart1[n]++; else numFoundOfThisPart2[n]++; } whichPart = n; } else { if (numTreePartsFound+1 > numPartsAllocated) { numPartsAllocated += 500; if (RealloateBits () == ERROR) return (ERROR); } x = &treePartsFound[numTreePartsFound * taxonLongsNeeded]; for (i=0; i 1) { for (i=0; i 0.0) { if (numFoundOfThisPart[whichPart] == 1) { aBrlens[whichPart] = bl; sBrlens[whichPart] = 0.0; } else { f = aBrlens[whichPart]; aBrlens[whichPart] += (bl - aBrlens[whichPart]) / (MrBFlt) (numFoundOfThisPart[whichPart]); sBrlens[whichPart] += (bl - aBrlens[whichPart]) * (bl - f); } if (sumtParams.numRuns > 1) { if (numFoundOfThisPart[whichPart] == 1) { aWithinBrlens[whichPart] = bl; sWithinBrlens[whichPart] = 0.0; sumB[whichPart] = bl; sumsqB[whichPart] = bl * bl; } else if (numFoundInRunOfPart[runIndex][whichPart] == 1) { aWithinBrlens[whichPart] = bl; sumB[whichPart] += bl; sumsqB[whichPart] += bl * bl; } else { f = aWithinBrlens[whichPart]; aWithinBrlens[whichPart] += (bl - aWithinBrlens[whichPart]) / (MrBFlt) (numFoundInRunOfPart[runIndex][whichPart]); sWithinBrlens[whichPart] += (bl - aWithinBrlens[whichPart]) * (bl - f); sumB[whichPart] += (aWithinBrlens[whichPart] - f); sumsqB[whichPart] += ((aWithinBrlens[whichPart] * aWithinBrlens[whichPart]) - (f * f)); } } } return (NO_ERROR); } int PrintBrlensToFile (void) { int i; /* print header */ for (i=0; ileft == NULL && p->right == NULL) { if (CheckString (p->label, taxaNames, &whichTaxon) == ERROR) { MrBayesPrint ("%s Could not find taxon %s in original list of taxa\n", spacer, p->label); goto errorExit; } whichTaxon--; if (taxaInfo[whichTaxon].isDeleted == YES) { prunedTaxa[whichTaxon] = YES; q = p->anc; if (q->left == p) sis = q->right; else sis = q->left; sis->length += q->length; if (q->anc == NULL) { MrBayesPrint ("%s Could not find root of tree\n", spacer); goto errorExit; } else qAnc = q->anc; if (qAnc->left == q) { qAnc->left = sis; sis->anc = qAnc; } else { qAnc->right = sis; sis->anc = qAnc; } p->left = p->right = p->anc = NULL; q->left = q->right = q->anc = NULL; numSumtTaxa--; deletedOne = YES; break; } } } if (deletedOne == YES) { i = 0; GetSumtDownPass (sumtRoot, downPass, &i); nNodes = i; } } while (deletedOne == YES); /* clean up on way out of function */ free (downPass); /* deroot tree, if we rooted it earlier */ if (wasDerooted == YES) { if (DerootSumtTree (sumtRoot, numSumtTaxa, outGroupNum) == ERROR) goto errorExit; } # if 0 if (isSumtTreeRooted == YES) MrBayesPrint ("%s Rooted tree %d:\n", spacer, numSumtTrees); else MrBayesPrint ("%s Unrooted tree %d:\n", spacer, numSumtTrees); ShowNodes (sumtRoot, 3, isSumtTreeRooted); # endif return (NO_ERROR); errorExit: if (downPass) free (downPass); return(ERROR); } int RealloateBits (void) { int i; numPartsAllocated += 100; if (memAllocs[ALLOC_TREEPARTS] == NO) goto errorExit; treePartsFound = (long *)realloc(treePartsFound, (size_t) (taxonLongsNeeded * numPartsAllocated * sizeof(long))); if (!treePartsFound) { MrBayesPrint ("%s Problem reallocating treePartsFound (%d)\n", spacer, taxonLongsNeeded * numPartsAllocated * sizeof(long)); goto errorExit; } if (memAllocs[ALLOC_NUMOFPART] == NO) goto errorExit; numFoundOfThisPart = (int *)realloc(numFoundOfThisPart, (size_t) (numPartsAllocated * sizeof(int))); if (!numFoundOfThisPart) { MrBayesPrint ("%s Problem reallocating numFoundOfThisPart (%d)\n", spacer, numPartsAllocated * sizeof(int)); goto errorExit; } if (sumtParams.numRuns > 1) { if (memAllocs[ALLOC_NUMINRUNOFPART] == NO) goto errorExit; for (i=0; i 1 && sumtBrlensDef == YES) { if (memAllocs[ALLOC_A_WITHIN_BRLENS] == NO) goto errorExit; aWithinBrlens = (MrBFlt *)realloc(aWithinBrlens, (size_t) (numPartsAllocated * sizeof(MrBFlt))); if (!aWithinBrlens) { MrBayesPrint ("%s Problem reallocating aWithinBrlens (%d bytes)\n", spacer, numPartsAllocated * sizeof(MrBFlt)); goto errorExit; } if (memAllocs[ALLOC_S_WITHIN_BRLENS] == NO) goto errorExit; sWithinBrlens = (MrBFlt *)realloc(sWithinBrlens, (size_t) (numPartsAllocated * sizeof(MrBFlt))); if (!sWithinBrlens) { MrBayesPrint ("%s Problem reallocating sWithinBrlens (%d bytes)\n", spacer, numPartsAllocated * sizeof(MrBFlt)); goto errorExit; } if (memAllocs[ALLOC_SUMB] == NO) goto errorExit; sumB = (MrBFlt *)realloc(sumB, (size_t) (numPartsAllocated * sizeof(MrBFlt))); if (!sumB) { MrBayesPrint ("%s Problem reallocating sumB (%d bytes)\n", spacer, numPartsAllocated * sizeof(MrBFlt)); goto errorExit; } if (memAllocs[ALLOC_SUMSQB] == NO) goto errorExit; sumsqB = (MrBFlt *)realloc(sumsqB, (size_t) (numPartsAllocated * sizeof(MrBFlt))); if (!sumsqB) { MrBayesPrint ("%s Problem reallocating sumsqB (%d bytes)\n", spacer, numPartsAllocated * sizeof(MrBFlt)); goto errorExit; } } return (NO_ERROR); errorExit: FreeBits(); return (ERROR); } int RealloateFullCompTrees (int whichList) { numFullCompTreesAllocated[whichList] += 100; if (memAllocs[ALLOC_FULLCOMPTREEINFO] == NO) goto errorExit; if (whichList == 0) { fullCompTreePartIds1 = (int *)realloc(fullCompTreePartIds1, (size_t) (2 * numTaxa * numFullCompTreesAllocated[whichList] * sizeof(int))); if (!fullCompTreePartIds1) { MrBayesPrint ("%s Problem reallocating fullCompTreePartIds1 (%d)\n", spacer, 2 * numTaxa * numFullCompTreesAllocated[whichList] * sizeof(int)); goto errorExit; } fullCompTreePartLengths1 = (MrBFlt *)realloc(fullCompTreePartLengths1, (size_t) (2 * numTaxa * numFullCompTreesAllocated[whichList] * sizeof(MrBFlt))); if (!fullCompTreePartLengths1) { MrBayesPrint ("%s Problem reallocating fullCompTreePartLengths1 (%d)\n", spacer, 2 * numTaxa * numFullCompTreesAllocated[whichList] * sizeof(MrBFlt)); goto errorExit; } } else { fullCompTreePartIds2 = (int *)realloc(fullCompTreePartIds2, (size_t) (2 * numTaxa * numFullCompTreesAllocated[whichList] * sizeof(int))); if (!fullCompTreePartIds2) { MrBayesPrint ("%s Problem reallocating fullCompTreePartIds2 (%d)\n", spacer, 2 * numTaxa * numFullCompTreesAllocated[whichList] * sizeof(int)); goto errorExit; } fullCompTreePartLengths2 = (MrBFlt *)realloc(fullCompTreePartLengths2, (size_t) (2 * numTaxa * numFullCompTreesAllocated[whichList] * sizeof(MrBFlt))); if (!fullCompTreePartLengths2) { MrBayesPrint ("%s Problem reallocating fullCompTreePartLengths2 (%d)\n", spacer, 2 * numTaxa * numFullCompTreesAllocated[whichList] * sizeof(MrBFlt)); goto errorExit; } } return (NO_ERROR); errorExit: FreeBits(); return (ERROR); } int RealloateFullTrees (void) { numFullTreesAllocated += 100; if (memAllocs[ALLOC_FULLTREEINFO] == NO) goto errorExit; numOfThisFullTree = (int *)realloc(numOfThisFullTree, (size_t) (numFullTreesAllocated * sizeof(int))); if (!numOfThisFullTree) { MrBayesPrint ("%s Problem reallocating numOfThisFullTree (%d)\n", spacer, numFullTreesAllocated * sizeof(int)); goto errorExit; } fullTreePartIds = (int *)realloc(fullTreePartIds, (size_t) (2 * numTaxa * numFullTreesAllocated * sizeof(int))); if (!fullTreePartIds) { MrBayesPrint ("%s Problem reallocating fullTreePartIds (%d)\n", spacer, 2 * numTaxa * numFullTreesAllocated * sizeof(int)); goto errorExit; } return (NO_ERROR); errorExit: FreeBits(); return (ERROR); } int ReorderParts (void) { int n, i, j; long *x, y, z, *newBits; newBits = (long *)SafeMalloc((size_t) (taxonLongsNeeded * sizeof(long))); if (!newBits) { MrBayesPrint ("%s Could not allocate newBits\n", spacer); return (ERROR); } numIncludedTaxa = 0; for (i=0; imemoryIndex] = YES; } j = 0; for (i=0; i<2*numTaxa; i++) { if (usedMemIndex[i] == NO) { if (j <= 1) { availableMemIndex[j] = i; j++; } else { /* some taxa are not included in the trees */ } } } first = &sumtNodes[availableMemIndex[0]]; second = &sumtNodes[availableMemIndex[1]]; /* root tree with previously down taxon as first right */ lft = sumtRoot->left; rht = sumtRoot; tempBrLen = lft->length; if (tempBrLen <= 0.0) tempBrLen = 0.0; lft->anc = rht->anc = first; rht->left = rht->right = NULL; first->left = lft; first->right = rht; first->anc = second; second->left = first; second->right = second->anc = NULL; lft->length = rht->length = tempBrLen * (MrBFlt) 0.5; first->length = second->length = 0.0; sumtRoot = second; isSumtTreeRooted = YES; /* update downpass sequence */ i = 0; GetSumtDownPass (sumtRoot, downPass, &i); nNodes = i; /* now, bring the outgroup around to the first right position */ isMarked = NO; for (i=0; ileft == NULL && p->right == NULL) { localTaxaFound[p->index] = YES; if (p->index == out) { p->marked = YES; isMarked = YES; } else p->marked = NO; } } if (isMarked == NO) { /* We will arbitrarily root it by the first taxon in the matrix that was actually found on the tree */ for (i=0; ileft == NULL && p->right == NULL) { localTaxaFound[p->index] = YES; if (p->index == sumtOut) { p->marked = YES; isMarked = YES; } else p->marked = NO; } } if (isMarked == NO) { MrBayesPrint ("%s Could not find outgroup taxon\n", spacer); goto errorExit; } localOutgroupNum = sumtOut; } for (i=0; ileft != NULL && p->right != NULL) if (p->left->marked == YES || p->right->marked == YES) p->marked = YES; } lft = sumtRoot->left; while (lft->left->index != localOutgroupNum && lft->right->index != localOutgroupNum) { if (lft->left->marked == YES && lft->right->marked == NO) { m1 = lft->left; um1 = lft->right; if (m1->left != NULL && m1->right != NULL) { if (m1->left->marked == YES) { m2 = m1->left; um2 = m1->right; } else { m2 = m1->right; um2 = m1->left; } lft->left = m2; lft->right = m1; m2->anc = m1->anc = lft; m1->left = um2; m1->right = um1; um1->anc = um2->anc = m1; m1->marked = NO; um1->length += m1->length; m2->length *= 0.5; m1->length = m2->length; } else { MrBayesPrint ("%s Rooting routine is lost (4)\n", spacer); goto errorExit; } } else if (lft->left->marked == NO && lft->right->marked == YES) { m1 = lft->right; um1 = lft->left; if (m1->left != NULL && m1->right != NULL) { if (m1->left->marked == YES) { m2 = m1->left; um2 = m1->right; } else { m2 = m1->right; um2 = m1->left; } lft->left = m1; lft->right = m2; m2->anc = m1->anc = lft; m1->left = um1; m1->right = um2; um1->anc = um2->anc = m1; m1->marked = NO; um1->length += m1->length; m2->length *= 0.5; m1->length = m2->length; } else { MrBayesPrint ("%s Rooting routine is lost (5)\n", spacer); goto errorExit; } } else { MrBayesPrint ("%s Rooting routine is lost (6)\n", spacer); goto errorExit; } } /* make certain outgroup is to the right of the root */ if (sumtRoot->left->left->index == localOutgroupNum) { m1 = sumtRoot->left->left; m2 = sumtRoot->left->right; lft = sumtRoot->left; lft->left = m2; lft->right = m1; } /* reindex internal nodes of tree */ i = numTaxa; FinishSumtTree (sumtRoot, &i, YES); /* free memory */ free (downPass); free (usedMemIndex); free (localTaxaFound); return (NO_ERROR); errorExit: if (downPass) free (downPass); if (usedMemIndex) free (usedMemIndex); if (localTaxaFound) free (localTaxaFound); return (ERROR); } void ShowBits (long *p, int nBitsToShow) { int i; long x, y; for (i=0; if = 0.0; for (i=nNodes-2; i>=0; i--) { p = allDownPass[i]; p->f = p->anc->f + p->length; if (p->left == NULL) { f = p->f / (screenWidth - strlen(p->label) - 2); if (f > scale) { scale = f; treeWidth = screenWidth - (int) strlen(p->label) - 2; } } } /* calculate x coordinates */ for (i=0; ix = (int) (0.5 + (p->f / scale)); } /* calculate y coordinates and lines to print */ for (i=nLines=0; ileft != NULL) { /* internal node */ for (q=p->left->sib; q->sib!=NULL; q=q->sib) ; p->y = (int) (0.5 + ((p->left->y + q->y) / 2.0)); } else { /* terminal node */ p->y = nLines; nLines += 2; } } /* print tree line by line */ for (i=0; iy != i) continue; /* this branch should be printed */ /* add branch */ if (p->anc == NULL) { /* this is the root of the whole tree */ printLine[p->x] = '+'; } else { /* this is an ordinary branch */ to = p->x; from = p->anc->x; for (k=from+1; k<=to; k++) printLine[k] = '-'; if (p == p->anc->left) { if (markLine[from] == 0) printLine[from] = '/'; else printLine[from] = '|'; markLine[from] ++; } else if (p->sib == NULL) { if (markLine[from] == 1) printLine[from] = '\\'; else printLine[from] = '|'; markLine[from] --; } if (p->left!=NULL) { if (from != to) printLine[to] = '+'; else printLine[to] = '|'; } else { /* add label if the branch is terminal */ sprintf(printLine+to+2,"%s",p->label); } } } /* check for cross branches */ for (j=0; j= 1 && printLine[j] == ' ') printLine[j] = '|'; } MrBayesPrintf (fp, "%s\n",printLine); } /* print scale */ k = (int) (floor (log10 (scale * 80))); scaleBar = pow (10, k); barLength = (int) (scaleBar / scale); if (barLength > 80) { barLength /= 10; scaleBar /= 10.0; } else if (barLength > 40) { barLength /= 5; scaleBar /= 5.0; } else if (barLength > 16) { barLength /= 2; scaleBar /= 2.0; } MrBayesPrintf (fp, "%s |", spacer); for (i=0; iindex); for (i=0; i<2*numTaxa; i++) { p = &conNodes[i]; if (!(p->left == NULL && p->sib == NULL && p->anc == NULL)) { printf ("%4d -- %2d ", i, p->index); if (p->sib != NULL) printf ("(%2d ", p->sib->index); else printf ("(-- "); if (p->left != NULL) printf ("%2d ", p->left->index); else printf ("-- "); if (p->anc != NULL) printf ("%2d)", p->anc->index); else printf ("--)"); if (p->left == NULL && p->anc != NULL) printf (" %s (%d)\n", p->label, p->x); else printf (" (%d)\n", p->x); } } return; } int ShowConTree (FILE *fp, int nNodes, PolyNode *root, int screenWidth, int showSupport) { int i, j, k, treeWidth, minBranchLength, maxWidth, isTreeDivided, printWidth, nLines, nodesToBePrinted, from, to, maxLabelLength; char *printLine, *markLine, temp[20]; PolyNode *p=NULL, *q, **allDownPass; minBranchLength = 5; isTreeDivided = NO; /* allocate space for printLine and markLine */ printLine = (char *) calloc ((2*screenWidth+1),sizeof(char)); if (!printLine) return ERROR; markLine = printLine + screenWidth + 1; /* allocate allDownPass */ allDownPass = (PolyNode **) SafeMalloc (nNodes * sizeof(PolyNode *)); if (!allDownPass) { free (printLine); return ERROR; } /* get fresh downpass sequence and internal node indices */ i = 0; GetConDownPass(allDownPass, root, &i); for (i=k=0; ileft == NULL) k++; for (i=0; ileft != NULL) p->index = k++; } /* calculate max length of labels initially set to hold largest node index number plus parentheses*/ maxLabelLength = (int) (3.0 + log10((MrBFlt)nNodes)); for (i=0; ileft == NULL && (int)strlen(p->label)>maxLabelLength) maxLabelLength = (int) strlen(p->label); } /* calculate remaining screen width for tree and maxWidth in terms of branches */ treeWidth = screenWidth - maxLabelLength - 1; maxWidth = treeWidth / minBranchLength; /* unmark whole tree */ for (i=0; imark = 0; nodesToBePrinted = nNodes; while (nodesToBePrinted > 0) { /* count depth of nodes in unprinted tree */ for (i=0; imark == 0) { p->x = 0; if (p->left != NULL && p->left->mark == 0) { for (q = p->left; q!=NULL; q=q->sib) { if (q->x > p->x) p->x = q->x; } p->x++; /* break when root of print subtree has been found */ if (p->x >= maxWidth) break; } } } /* if internal node then find largest nonprinted subtree among descendant nodes */ if (p->anc != NULL) { for (q=p->left; q!=NULL; q=q->sib) { if (q->x == p->x - 1 && q->mark == 0) p = q; } MrBayesPrintf (fp, "%s Subtree rooted at node %d:\n\n", spacer, p->index); isTreeDivided = YES; } else if (isTreeDivided == YES) MrBayesPrintf (fp, "%s Root part of tree:\n\n", spacer); /* mark subtree for printing and translate x coordinates from depth to position */ if (p->anc == NULL) printWidth = p->x; else printWidth = p->x + 1; p->mark = 1; p->x = (int) (treeWidth - 0.5 - ((treeWidth - 1) * (p->x / (MrBFlt) printWidth))); for (i=nNodes-2; i>=0; i--) { p = allDownPass[i]; if (p->mark == 0 && p->anc->mark == 1) { p->mark = 1; p->x = (int) (treeWidth - 0.5 - ((treeWidth - 1) * (p->x / (MrBFlt) printWidth))); } } /* calculate y coordinates of nodes to be printed and lines to print */ for (i=nLines=0; imark == 1) { if (p->left != NULL && p->left->mark == 1) { /* internal node */ for (q=p->left->sib; q->sib!=NULL; q=q->sib) ; p->y = (int) (0.5 + ((p->left->y + q->y) / 2.0)); } else { /* terminal node */ p->y = nLines; nLines += 2; } } } /* print subtree line by line */ for (i=0; imark != 1 || p->y != i) continue; /* this branch should be printed add label if the branch is terminal in tree to be printed */ if (p->left == NULL) sprintf(printLine+treeWidth+1,"%s", p->label); else if (p->left->mark == 2) sprintf(printLine+treeWidth+1,"(%d)", p->index); /* add branch */ if (p->anc == NULL) { /* this is the root of the whole tree */ printLine[p->x] = '+'; nodesToBePrinted--; } else if (p->anc->mark == 0) { /* this is a root of a subtree this branch will have to be printed again so do not decrease nodesToBePrinted */ to = p->x; from = 0; for (k=from; ksupport + 0.5)); else *temp='\0'; from = (int)(from + 1.5 + ((to - from - 1 - strlen(temp)) / 2.0)); for (k=0; temp[k]!='\0'; k++) printLine[from++] = temp[k]; } else { /* this is an ordinary branch */ to = p->x; from = p->anc->x; for (k=from+1; k<=to; k++) printLine[k] = '-'; if (p == p->anc->left) { printLine[from] = '/'; markLine[from] = 1; } else if (p->sib == NULL) { printLine[from] = '\\'; markLine[from] = 0; } if (p->left!=NULL && p->left->mark!=2) { printLine[to] = '+'; if (showSupport == YES) sprintf(temp, "%d", (int) (p->support + 0.5)); else *temp='\0'; from = (int)(from + 1.5 + ((to - from - 1 - strlen(temp)) / 2.0)); for (k=0; temp[k]!='\0'; k++) printLine[from++] = temp[k]; } nodesToBePrinted--; } } /* check for cross branches */ for (j=0; jmark == 1) { if (p->anc == NULL) p->mark = 2; else if (p->anc->mark == 0) p->mark = 0; /* this branch will have to be printed again */ else p->mark = 2; } } } /* next subtree */ free (allDownPass); free (printLine); return NO_ERROR; } void ShowParts (long *p, int nTaxaToShow) { int i, flipBits; long x, y; flipBits = NO; if (isSumtTreeRooted == NO) { y = p[0]; x = 1 << (0 % nBitsInALong); if (x & y) flipBits = YES; } for (i=0; ileft == NULL && p->right == NULL && p->anc != NULL) { MrBayesPrint("%*cN %d (l=%d r=%d a=%d) %lf (%s) ", indent, ' ', SumtDex(p), SumtDex(p->left), SumtDex(p->right), SumtDex(p->anc), p->length, p->label); } else if (p->left != NULL && p->right == NULL && p->anc == NULL) { if (isThisTreeRooted == NO) { if (p->label[0] == '\0' || p->label[0] == '\n' || p->label[0] == ' ') MrBayesPrint("%*cN %d (l=%d r=%d a=%d) (---) ", indent, ' ', SumtDex(p), SumtDex(p->left), SumtDex(p->right), SumtDex(p->anc)); else MrBayesPrint("%*cN %d (l=%d r=%d a=%d) (%s) ", indent, ' ', SumtDex(p), SumtDex(p->left), SumtDex(p->right), SumtDex(p->anc), p->label); } else { MrBayesPrint("%*cN %d (l=%d r=%d a=%d) X.XXXXXX ", indent, ' ', SumtDex(p), SumtDex(p->left), SumtDex(p->right), SumtDex(p->anc)); } } else { if (p->anc != NULL) { if (p->anc->anc == NULL && isThisTreeRooted == YES) MrBayesPrint("%*cN %d (l=%d r=%d a=%d) X.XXXXXX ", indent, ' ', SumtDex(p), SumtDex(p->left), SumtDex(p->right), SumtDex(p->anc)); else MrBayesPrint("%*cN %d (l=%d r=%d a=%d) %lf ", indent, ' ', SumtDex(p), SumtDex(p->left), SumtDex(p->right), SumtDex(p->anc), p->length); } } MrBayesPrint ("\n"); ShowSumtNodes (p->left, indent + 2, isThisTreeRooted); ShowSumtNodes (p->right, indent + 2, isThisTreeRooted); } } void SortInts (int *item, int *assoc, int count, int descendingOrder) { SortInts2 (item, assoc, 0, count-1, descendingOrder); } void SortInts2 (int *item, int *assoc, int left, int right, int descendingOrder) { register int i, j, x, y; if (descendingOrder == YES) { i = left; j = right; x = item[(left+right)/2]; do { while (item[i] > x && i < right) i++; while (x > item[j] && j > left) j--; if (i <= j) { y = item[i]; item[i] = item[j]; item[j] = y; if (assoc) { y = assoc[i]; assoc[i] = assoc[j]; assoc[j] = y; } i++; j--; } } while (i <= j); if (left < j) SortInts2 (item, assoc, left, j, descendingOrder); if (i < right) SortInts2 (item, assoc, i, right, descendingOrder); } else { i = left; j = right; x = item[(left+right)/2]; do { while (item[i] < x && i < right) i++; while (x < item[j] && j > left) j--; if (i <= j) { y = item[i]; item[i] = item[j]; item[j] = y; if (assoc) { y = assoc[i]; assoc[i] = assoc[j]; assoc[j] = y; } i++; j--; } } while (i <= j); if (left < j) SortInts2 (item, assoc, left, j, descendingOrder); if (i < right) SortInts2 (item, assoc, i, right, descendingOrder); } } void SortIndParts (int *item, MrBFlt *assoc, int count, int descendingOrder) { SortIndParts2 (item, assoc, 0, count-1, descendingOrder); } void SortIndParts2 (int *item, MrBFlt *assoc, int left, int right, int descendingOrder) { register int i, j, x; MrBFlt y; if (descendingOrder == YES) { i = left; j = right; x = item[(left+right)/2]; do { while (item[i] > x && i < right) i++; while (x > item[j] && j > left) j--; if (i <= j) { y = (MrBFlt) item[i]; item[i] = item[j]; item[j] = (int) y; if (assoc) { y = assoc[i]; assoc[i] = assoc[j]; assoc[j] = y; } i++; j--; } } while (i <= j); if (left < j) SortIndParts2 (item, assoc, left, j, descendingOrder); if (i < right) SortIndParts2 (item, assoc, i, right, descendingOrder); } else { i = left; j = right; x = item[(left+right)/2]; do { while (item[i] < x && i < right) i++; while (x < item[j] && j > left) j--; if (i <= j) { y = (MrBFlt) item[i]; item[i] = item[j]; item[j] = (int) y; if (assoc) { y = assoc[i]; assoc[i] = assoc[j]; assoc[j] = y; } i++; j--; } } while (i <= j); if (left < j) SortIndParts2 (item, assoc, left, j, descendingOrder); if (i < right) SortIndParts2 (item, assoc, i, right, descendingOrder); } } int SortParts (int *item, int count) { int i, *tempVect; SortParts2 (item, 0, count-1); tempVect = (int *)SafeMalloc((size_t) (numTreePartsFound * sizeof(int))); if (!tempVect) { MrBayesPrint ("%s Problem allocating tempVect (%d)\n", spacer, numTreePartsFound * sizeof(int)); goto errorExit; } for (i=0; i %d\n", partOrigOrder[i], i);*/ free (tempVect); return (NO_ERROR); errorExit: if (tempVect) free (tempVect); return (ERROR); } void SortParts2 (int *item, int left, int right) { register int i, j, k; int yI; long yL; MrBFlt x, y; i = left; j = right; x = (MrBFlt) item[(left+right)/2]; do { while (item[i] > x && i < right) i++; while (x > item[j] && j > left) j--; if (i <= j) { yI = item[i]; item[i] = item[j]; item[j] = yI; yI = partOrigOrder[i]; partOrigOrder[i] = partOrigOrder[j]; partOrigOrder[j] = yI; if (sumtBrlensDef == YES) { y = aBrlens[i]; aBrlens[i] = aBrlens[j]; aBrlens[j] = y; y = sBrlens[i]; sBrlens[i] = sBrlens[j]; sBrlens[j] = y; } if (sumtParams.numRuns > 1 && sumtBrlensDef == YES) { y = aWithinBrlens[i]; aWithinBrlens[i] = aWithinBrlens[j]; aWithinBrlens[j] = y; y = sWithinBrlens[i]; sWithinBrlens[i] = sWithinBrlens[j]; sWithinBrlens[j] = y; y = sumB[i]; sumB[i] = sumB[j]; sumB[j] = y; y = sumsqB[i]; sumsqB[i] = sumsqB[j]; sumsqB[j] = y; } for (k = 0; k 1) { for (k=0; kindex; } int TreeProb (void) { int i, j, n, num, targetNode, nBits, nextConNode, isCompat, localOutgroupNum, origPartNum, reorderedPartNum, *tempTreeNum=NULL, *tempNumOfTree=NULL, nInSets[5]; long x, *mask, *partition, *ingroupPartition, *outgroupPartition=NULL; MrBFlt treeProb, cumTreeProb; char tempName[100]; PolyNode *cp, *q, *r, *ql, *pl; /* check if we need to do this */ if (sumtParams.calcTrprobs == NO) return (NO_ERROR); MrBayesPrint ("%s Calculating tree probabilities...\n\n", spacer); /* check that we have at least three species */ j = 0; for (i=0; ix counts number of subtended terminals make sure conRoot->left is in outgroup */ conRoot = &conNodes[numIncludedTaxa]; conRoot->anc = conRoot->sib = NULL; conRoot->x = numIncludedTaxa; j = FirstTaxonInPartition(outgroupPartition, taxonLongsNeeded); conRoot->left = cp = &conNodes[j]; cp->anc = conRoot; cp->x = 1; for (i=0; isib = &conNodes[i]; cp = cp->sib; cp->anc = conRoot; cp->x = 1; } } cp->sib = NULL; /* Resolve bush according to partitions. Partitions may include incompatible ones. */ nextConNode = numIncludedTaxa + 1; if (isSumtTreeRooted == YES) targetNode = 2 * numIncludedTaxa - 2; else targetNode = 2 * numIncludedTaxa - 3; for (i=0; i 1) /* this is an informative partition */ { /* find anc of partition */ j = FirstTaxonInPartition (partition, taxonLongsNeeded); for (cp = &conNodes[j]; cp!=NULL; cp = cp->anc) if (cp->x > nBits) break; /* Do not include if incompatible with ancestor or any of descendants. Do not check terminals or root because it is redundant and partitions have not been set for those. */ isCompat = YES; if (cp->anc != NULL && !IsPartCompatible(partition, cp->partition, taxonLongsNeeded)) isCompat = NO; for (q=cp->left; q!=NULL; q=q->sib) { if (q->x > 1 && !IsPartCompatible(q->partition, partition, taxonLongsNeeded)) isCompat = NO; if (isCompat == NO) { MrBayesPrint ("%s Found an incompatible partition in the tree (1 %d %d)\n", spacer, num, n); goto errorExit; } } if (isCompat == NO) { MrBayesPrint ("%s Found an incompatible partition in the tree (2)\n", spacer); goto errorExit; } /* check for number of nodes in tree */ if (nextConNode > targetNode) { MrBayesPrint ("%s Too many nodes in tree\n", spacer); goto errorExit; } /* set new node */ q = &conNodes[nextConNode++]; q->x = nBits; q->partition = partition; /* go through descendants of anc */ ql = pl = NULL; for (r=cp->left; r!=NULL; r=r ->sib) { /* test if r is in the new partition or not */ if ((r->x > 1 && IsPartNested(r->partition, partition, taxonLongsNeeded)) || (r->x == 1 && (partition[r->index / nBitsInALong] & (1 << (r->index % nBitsInALong))) != 0)) { /* r is in the partition */ if (ql == NULL) q->left = r; else ql->sib = r; ql = r; r->anc = q; } else { /* r is not in the partition */ if (pl == NULL) cp->left = r; else pl->sib = r; pl = r; } } /* terminate new sib-node chain */ ql->sib = NULL; /* new node is last in old sib-node chain */ pl->sib = q; q->sib = NULL; q->anc = cp; } else /* singleton partition */ { j = FirstTaxonInPartition(partition, taxonLongsNeeded); q = &conNodes[j]; } } /* get probability of tree */ treeProb = (MrBFlt)tempNumOfTree[num] / (MrBFlt) (sumtParams.numRuns * numSumTreesSampled); cumTreeProb += treeProb; if (cumTreeProb >= 0.0 && cumTreeProb < 0.5) nInSets[0]++; else if (cumTreeProb >= 0.5 && cumTreeProb < 0.9) nInSets[1]++; else if (cumTreeProb >= 0.9 && cumTreeProb < 0.95) nInSets[2]++; else if (cumTreeProb >= 0.95 && cumTreeProb < 0.99) nInSets[3]++; else nInSets[4]++; /* draw tree to stdout */ if (sumtParams.showSumtTrees == YES) { MrBayesPrint ("\n%s Tree %d (p = %1.3lf, P = %1.3lf):\n\n", spacer, num+1, treeProb, cumTreeProb); ShowConTree (stdout, nextConNode, conRoot, 80, NO); } /* draw tree to file */ if (num == 0) { MrBayesPrintf (fpTrees, "[This file contains the trees that were found during the MCMC\n"); MrBayesPrintf (fpTrees, "search, sorted by posterior probability. \"p\" indicates the\n"); MrBayesPrintf (fpTrees, "posterior probability of the tree whereas \"P\" indicates the\n"); MrBayesPrintf (fpTrees, "cumulative posterior probability.]\n\n"); MrBayesPrintf (fpTrees, "begin trees;\n"); MrBayesPrintf (fpTrees, " translate\n"); j = 0; for (i=0; ianc != NULL) if (p->anc->left == p) fprintf (fp, "("); for (q = p->left; q != NULL; q = q->sib) { if (q->anc->left != q) /* Note that q->anc always exists (it is p) */ fprintf (fp, ","); WriteConTree (q, fp, showSupport); } if (p->left == NULL) { if (sumtBrlensDef == YES) fprintf (fp, "%s:%lf", p->label, p->length); else fprintf (fp, "%s", p->label); } if (p->sib == NULL && p->anc != NULL) { if (p->anc->anc != NULL) { if (sumtBrlensDef == YES && showSupport == NO) fprintf (fp, "):%lf", p->anc->length); else if (sumtBrlensDef == NO && showSupport == YES) fprintf (fp, ")%1.2lf", p->anc->support/100.0); else if (sumtBrlensDef == YES && showSupport == YES) fprintf (fp, ")%1.2lf:%lf", p->anc->support/100.0, p->anc->length); else fprintf (fp, ")"); } else fprintf (fp, ")"); } } void WriteTree (PolyNode *p, FILE *fp) { PolyNode *q; if (p->anc != NULL) if (p->anc->left == p) fprintf (fp, "("); for (q = p->left; q != NULL; q = q->sib) { if (q->anc->left != q) /* Note that q->anc always exists (it is p) */ fprintf (fp, ","); WriteTree (q, fp); } if (p->left == NULL) { fprintf (fp, "%d", p->index+1); } if (p->sib == NULL && p->anc != NULL) { fprintf (fp, ")"); } }