/* * MrBayes 2.0 * * John P. Huelsenbeck * Department of Biology * University of Rochester * Rochester, NY 14627 * * johnh@brahms.biology.rochester.edu * * Fredrik Ronquist * Dept. Systematic Zoology * Evolutionary Biology Centre * Uppsala University * Norbyv. 18D, SE-752 36 Uppsala, Sweden * * fredrik.ronquist@ebc.uu.se * */ #include #include #include #include #include #include #include "jph.h" #include "tree.h" #include "command.h" #include "mcmc.h" #define NUMCOMMANDS 176 #define NUMTOKENCHARS 14 #define HIERARCHYSIZE 610 #define SCREENWIDTH 60 #define SCREENWIDTH2 61 #define SCREENHEIGTH 15 #define HEADER_SIZE 1000 #define MAX_PARTITIONS 10000 #define NO_ACTION -1 #define ERROR_ACTION -2 #define EQUALSIGN 1 #define COLON 2 #define SEMICOLON 3 #define COMMA 4 #define POUNDSIGN 5 #define QUESTIONMARK 6 #define DASH 7 #define LEFTPAR 8 #define RIGHTPAR 9 #define LEFTCOMMENT 10 #define RIGHTCOMMENT 11 #define ALPHA 12 #define NUMBER 13 #define RETURNSYMBOL 14 #define ASTERISK 15 #define BACKSLASH 16 #define FORWARDSLASH 17 #define EXCLAMATIONMARK 18 #define PERCENT 19 #define AMPERSTAND 20 #define TILDE 21 #define PLUS 22 #define CARAT 23 #define DOLLARSIGN 24 #define ATSIGN 25 #define PIPESIGN 26 #define LEFTCURLY 27 #define RIGHTCURLY 28 #define UNKNOWN_TOKEN_TYPE 29 #define ANYTOKEN 30 #define CMD_EXECUTE 0 #define CMD_QUIT 1 #define CMD_ENTER_INPUT_FILE_NAME1 2 #define CMD_ENTER_INPUT_FILE_NAME2 3 #define CMD_INIT_NTAX 4 #define CMD_INIT_NCHAR 5 #define CMD_INIT_DATA_BLOCK 6 #define CMD_ALLOC_MATRIX 7 #define CMD_INIT_GAPCHAR 8 #define CMD_INIT_MATCHCHAR 9 #define CMD_INIT_MISSINGCHAR 10 #define CMD_INIT_INTERLEAVE 11 #define CMD_INIT_DATATYPE 12 #define CMD_EXIT_DATA_BLOCK 13 #define CMD_INIT_READ_MATRIX 14 #define CMD_FILL_MATRIX 15 #define CMD_FINISH_READ_MATRIX 16 #define CMD_ENTER_HELP_Q 17 #define CMD_GET_USER_SOME_HELP 18 #define CMD_INIT_HELP 19 #define CMD_SET_NST 20 #define CMD_SET_RATES 21 #define CMD_SET_BASE_FREQ 22 #define CMD_SET_CLOCK 23 #define CMD_END_NST 24 #define CMD_INIT_TRATIO 25 #define CMD_INIT_REVRATES 26 #define CMD_INIT_REVRATES2 27 #define CMD_INIT_REVRATES3 28 #define CMD_INIT_NONREVRATES 29 #define CMD_INIT_NONREVRATES2 30 #define CMD_INIT_NONREVRATES3 31 #define CMD_INIT_SHAPE 32 #define CMD_INIT_SITERATES 33 #define CMD_INIT_SITERATES2 34 #define CMD_INIT_PARTITION 35 #define CMD_INIT_CHARSET_NAME 36 #define CMD_INIT_CHARSET_NAME2 37 #define CMD_INIT_CHARSET 38 #define CMD_INIT_CHARSET2 39 #define CMD_ENTER_MRBAYES_BLOCK 40 #define CMD_EXIT_MRBAYES_BLOCK 41 #define CMD_INIT_EXCLUDE 42 #define CMD_EXCLUDE 43 #define CMD_EXCLUDE2 44 #define CMD_EXCLUDE3 45 #define CMD_EXCLUDE4 46 #define CMD_DEF_PARTITION 47 #define CMD_PARTITION 48 #define CMD_PARTITION2 49 #define CMD_PARTITION3 50 #define CMD_PARTITION4 51 #define CMD_PARTITION5 52 #define CMD_PARTITION6 53 #define CMD_DEF_OUTGROUP1 54 #define CMD_DEF_OUTGROUP2 55 #define CMD_SHOW_CHAR_STAT 56 #define CMD_CALIBRATION1 57 #define CMD_CALIBRATION2 58 #define CMD_CALIBRATION3 59 #define CMD_CALIBRATION4 60 #define CMD_CALIBRATION5 61 #define CMD_CALIBRATION6 62 #define CMD_CALIBRATION7 63 #define CMD_CONSTRAINT1 64 #define CMD_CONSTRAINT2 65 #define CMD_CONSTRAINT3 66 #define CMD_CONSTRAINT4 67 #define CMD_CONSTRAINT5 68 #define CMD_INIT_INCLUDE 69 #define CMD_INCLUDE 70 #define CMD_INCLUDE2 71 #define CMD_INCLUDE3 72 #define CMD_INCLUDE4 73 #define CMD_DELETE1 74 #define CMD_DELETE2 75 #define CMD_DELETE3 76 #define CMD_RESTORE1 77 #define CMD_RESTORE2 78 #define CMD_RESTORE3 79 #define CMD_MCMC 80 #define CMD_MCMC1 81 #define CMD_MCMC2 82 #define CMD_MCMC3 83 #define CMD_MCMC4 84 #define CMD_MCMC5 85 #define CMD_MCMC6 86 #define CMD_MCMC7 87 #define CMD_MCMC8 88 #define CMD_MCMC9 89 #define CMD_PRIOR1 90 #define CMD_PRIOR2 91 #define CMD_PRIOR3 92 #define CMD_PRIOR4 93 #define CMD_PRIOR5 94 #define CMD_PRIOR6 95 #define CMD_PRIOR7 96 #define CMD_PRIOR8 97 #define CMD_PRIOR9 98 #define CMD_PRIOR10 99 #define CMD_PRIOR11 100 #define CMD_PRIOR12 101 #define CMD_PRIOR13 102 #define CMD_PRIOR14 103 #define CMD_PRIOR15 104 #define CMD_PRIOR16 105 #define CMD_PRIOR17 106 #define CMD_PRIOR18 107 #define CMD_PRIOR19 108 #define CMD_PRIOR20 109 #define CMD_PRIOR21 110 #define CMD_PRIOR22 111 #define CMD_PRIOR23 112 #define CMD_PRIOR24 113 #define CMD_PRIOR25 114 #define CMD_PRIOR26 115 #define CMD_PRIOR27 116 #define CMD_TREE 117 #define CMD_TREE1 118 #define CMD_TREE2 119 #define CMD_TREE3 120 #define CMD_TREE4 121 #define CMD_TREE5 122 #define CMD_TREE6 123 #define CMD_TREE7 124 #define CMD_TREE8 125 #define CMD_SHOWUSERTREE 126 #define CMD_INIT_ENFORCECAL 127 #define CMD_INIT_ENFORCECON 128 #define CMD_INIT_INFERANC 129 #define CMD_INIT_NCAT 130 #define CMD_ROOT_TREE 131 #define CMD_DEROOT_TREE 132 #define CMD_SHOWNODES 133 #define CMD_SAVEBRLENS 134 #define CMD_SET_SEQ_ERROR 135 #define CMD_SET_SP_RATE 136 #define CMD_SET_EX_RATE 137 #define CMD_SET_SAMP_FRAC 138 #define CMD_SUMT_FILENAME 139 #define CMD_SUMT_BURN 140 #define CMD_SUMT_READFILE 141 #define CMD_TRANS_INIT 142 #define CMD_TRANS_NEW_NM 143 #define CMD_TRANS_OLD_NM 144 #define CMD_TRANS_COMMA 145 #define CMD_TRANS_END 146 #define CMD_SUM_TREE 147 #define CMD_SUM_TREE1 148 #define CMD_SUM_TREE2 149 #define CMD_SUM_TREE3 150 #define CMD_SUM_TREE4 151 #define CMD_SUM_TREE5 152 #define CMD_SUM_TREE6 153 #define CMD_SUM_TREE7 154 #define CMD_SUM_TREE8 155 #define CMD_INIT_TREE_BLK 156 #define CMD_SUM_TREE_END 157 #define CMD_SUMT_HALFCOMP 158 #define CMD_SUMT_ALLCOMP 159 #define CMD_SUMP_FILENAME 160 #define CMD_SUMP_BURN 161 #define CMD_SUMP_READFILE 162 #define CMD_SUM_PARM1 163 #define CMD_SUM_PARM2 164 #define CMD_SUM_PARM3 165 #define CMD_SUM_PARM4 166 #define CMD_SUM_PARM5 167 #define CMD_SUM_PARM6 168 #define CMD_SUM_PARM7 169 #define CMD_SUM_PARM8 170 #define CMD_SUM_PARM9 171 #define CMD_SUM_PARM10 172 #define CMD_AUTOCLOSE1 173 #define CMD_AUTOCLOSE2 174 #define CMD_NPERB 175 #define CMD_MCMCP 176 #define CMD_SET_AAMODEL1 177 #define CMD_SET_AAMODEL2 178 #define CMD_SET_AAMODEL3 179 #define CMD_SET_AAMODEL4 180 #define CMD_SET_AAMODEL5 181 #define CMD_SET_ANCFILE1 182 #define CMD_SET_ANCFILE2 183 #define CMD_SET_CODON1 184 #define CMD_SET_CODON2 185 #define CMD_SET_OMEGAVAR1 186 #define CMD_SET_OMEGAVAR2 187 #define CMD_SET_OMEGAVAR3 188 #define CMD_SET_OMEGAVAR4 189 #define CMD_SET_OMEGAVAR5 190 #define CMD_SET_OMEGA 191 #define CMD_OMEGAPR1 192 #define CMD_OMEGAPR2 193 #define CMD_OMEGAPR3 194 #define CMD_OMEGAPR4 195 #define CMD_OMEGAPR5 196 #define CMD_OMEGAPR6 197 #define CMD_SET_CODE 198 #define CMD_INFERSEL_YES 199 #define CMD_INFERSEL_NO 200 #define CMD_INIT_LAG1 201 #define CMD_INIT_LAG2 202 #define CMD_SET_LAG_PR1 203 #define CMD_SET_LAG_PR2 204 #define CMD_SET_LAG_PR3 205 #define CMD_SET_LAG_PR4 206 #define CMD_SET_INFERRATES1 207 #define CMD_SET_INFERRATES2 208 #define CMD_SET_COVARION1 209 #define CMD_SET_COVARION2 210 #define CMD_SET_PROPS1 211 #define CMD_SET_PROPS2 212 #define CMD_SET_COVTYPE 213 #define CMD_SET_CODING 214 #define CMD_SET_STATE_PR1 215 #define CMD_SET_STATE_PR2 216 #define CMD_SET_MOVE_PROPS 217 #define CMD_SET_MCMC_BURN 218 #define CMD_SET_CYCLEFREQ 219 #define CMD_SET_CALCPSEUDOBF 220 #define CMD_SET_PARS_CRIT 221 #define CMD_SET_TYPESET1 222 #define CMD_SET_TYPESET2 223 #define CMD_SET_TYPESET3 224 #define CMD_SET_TYPESET4 225 #define CMD_SET_TYPESET5 226 #define CMD_SET_TYPESET6 227 #define CMD_SET_TYPESET7 228 #define CMD_SET_TYPESET8 229 #define CMD_SET_TYPESET9 230 #define CMD_SET_TYPESET10 231 #define CMD_SET_TYPESET11 232 #define CMD_SET_STATESSET1 233 #define CMD_SET_STATESSET2 234 #define CMD_SET_STATESSET3 235 #define CMD_SET_STATESSET4 236 #define CMD_SET_STATESSET5 237 #define CMD_SET_STATESSET6 238 #define CMD_SET_STATESSET7 239 #define CMD_SET_STATESSET8 240 #define CMD_SET_STATESSET9 241 #define CMD_SET_STATESSET10 242 #define CMD_SET_STATESSET11 243 #define CMD_SET_ASSUME1 244 #define CMD_SET_ASSUME2 245 #define CONTREE typedef struct cmdNode { struct cmdNode *parent, *left, *sib, *rNde; int index, cmd, action; char cmdName[100]; } command_list; typedef struct conNode { struct conNode *left, *leftSib, *rightSib, *anc; int index; char label[100]; double length; } ConTreeNode; #if defined (CONTREE) typedef struct pnode { struct pnode *left, *sib, *anc; int x, y, mark, index; double length, support, f; char *label; unsigned int *partition; } PolyNode; #endif char ChangeCase (char c); void CleanExit (void); int DerootSumTree (TreeNode *r); int Dex(command_list *p); int Dex3 (ConTreeNode *p); void DisplayBits (unsigned value); int DoAction (int actionType); void ErrorExit (void); void GetToken (int *tokenType); void GetUserDownPass (TreeNode *p, TreeNode **x, int *y); void GetUserHelp (void); int IsCommandValid (int commdandNumber); int IsIn (char ch, char s[]); int IsWhite (char c); int NucID (char nuc); int ParseFile (void); int ParseToken (int tokenType); void PrepareToken (void); int ProtID (char aa); void ReadInTempToken (char tempToken[NUMTOKENCHARS], int i); int ReadNexusFile (void); int RealloateTreeParts (void); int ResID (char nuc); int RootSumTree (TreeNode *p); int SetBits (int n); int ShowTree (TreeNode *r, int isThisTreeRooted); void SortParts (int *item, int count); void SortParts2 (int *item, int left, int right); int StandID (char nuc); int SummarizeParams (void); int SummarizeTrees (void); # if defined (CONTREE) int FirstTaxonInPartition (unsigned int *partition, int length); void FlipBits(unsigned int *partition, int length, unsigned int lastMask); void GetConDownPass (PolyNode **downPass, PolyNode *p, int *i); int IsPartCompatible (unsigned int *smaller, unsigned int *larger, int length); int IsPartNested (unsigned int *smaller, unsigned int *larger, int length); int ShowConTree (FILE *fp, int nNodes, PolyNode *root, int screenWidth); int ShowConPhylogram (FILE *fp, int nNodes, PolyNode *root, int screenWidth); void WriteConTree (PolyNode *p, FILE *fp); # else void SortCompatParts (int *item, unsigned *item2, double *item3, int count); void SortCompatParts2 (int *item, unsigned *item2, double *item3, int left, int right); void WriteConTreeToScreen (ConTreeNode *p, FILE *fp); # endif extern int nTaxa, nChar, *origDataMatrix, isNtaxaInitialized, isNcharInitialized, isMatrixAllocated, dataType, *excludedSites, *sitePartitions, *charType, *charNStates, molecularClock, *calibrations, *constraints, nConstraints, nCalibrations, *excludedTaxa, nGen, sampleFreq, printFreq, numChains, isUserTreeDefined, isTreeRooted, enforceConstraints, enforceCalibrations, inferAncStates, nRateCats, isPartitionDefined, nPartitions, outgroupNum, saveBrlens, burnIn, consensusType, paramBurn, autoclose, nPerturbations, userBrlensDef, aaModelType, saveAncFiles, enforceCodonModel, aaCode, inferPoseSel, lagNumber, inferRates, useCovarion, mcmcBurn, cycleFreq, calcPseudoBf, useParsCriterion; extern long int randomSeed; extern char *taxonLabels, subModel[20], baseFreqModel[20], ratesModel[20], clockModel[20], partitionName[50], defPartitionName[50], tRatioModel[20], revRatesModel[20], nonRevRatesModel[20], shapeModel[20], siteRatesModel[20], calibrationLabels[32][100], constraintLabels[32][100], outFile[100], qmatprModel[20], basefreqprModel[20], brlenprModel[20], shapeprModel[20], siterateprModel[20], startingTreeModel[20], sumTreeFile[100], sumParmFile[100], omegaRatioModel[20], codonModelType[20], omegaprModel[20], lagModel[20], lagprModel[20], covarionModelType[20], codingType[20], stateFreqPrModel[20]; extern double revRates[6], nonRevRates[12], shape, tRatio, siteRates[50], calibrationTimes[32][2], chainTemp, qmatprExp, qmatprUni[2], basefreqprDir, brlenprExp, brlenprUni[2], shapeprExp, shapeprUni[2], siterateprExp, siterateprUni[2], seqErrorProb, speciationRate, extinctionRate, samplingFraction, omegaRatio, omegaprUni[2], omegaprExp, lagprUni[2], lagprExp, qMatWinProp, piAlphaProp, rateWinProp, tuningProp, clockTuneProp, timeTuneProp, mTuneProp, extendPProp, wormTuneProp, bdWinProp, omegaWinProp, rhoWinProp, omegaAlphaProp, invSlideProp, switchWinProp, alphaWinProp, nodeTimeTuneProp, reconLimProp, tbrExtendPProp, tbrTuneProp, blTuneProp, statefreqprDir, moveRates[100]; extern TreeNode *userTree, *userTreeRoot; static int expectedCommands[HIERARCHYSIZE], isInterleaved, *matrixCounter, returnEncountered, whichTaxon, isFirstInterleavedBlock, lineNumber, generalHelp, globalInt, *charSet, nCharSets, charSetFrom, charSetTo, charSetWait, charSetFreq, charSetSlashWait, charSetPotWait, whichPartition, isFirstNode, nextAvailableNode, intNodeNum, transNumber, isSumTreeRooted, isTranslated, numSumTrees, nUnsignedNeeded, numTreePartsFound, *numFoundOfThisPart, numSumTreesSampled, foundBrlens, nColumns, headerSize, tempNumColumns, curNumRows, whichRow, whichColumn, numPartsAllocated, changeToLower, *typeDef, nTypeDefs, charOrdering, *nStatesDef, nNStatesDefs, nCharStates; static unsigned *treeBits, *treePartsFound, *taxonMask, propSet; static double *aBrlens, *aOldBrlens, *sBrlens, *paramInfo; static command_list cmdList[HIERARCHYSIZE], *cmdRoot, *cmdPtr; static TreeNode *pPtr, *qPtr, *sumTree, *sumTreeRoot, **sumTreeDP; #if !defined (CONTREE) static ConTreeNode *constraintTree; #endif static char *tokenP, token[STRING_LENGTH], tokenType, inFile[100], *headerNames, gapChar, matchChar, missingChar, helpRequest[100], charSetLabels[50][100], typeDefLabels[50][100], nStatesDefLabels[50][100], *transNewNames, *transOldNames; static char validTokens[NUMCOMMANDS][NUMTOKENCHARS] = { "execute", /* 0 */ "begin", /* 1 */ "end", /* 2 */ "mrbayes", /* 3 */ "data", /* 4 */ "trees", /* 5 */ "help", /* 6 */ "quit", /* 7 */ "lset", /* 8 */ "nst", /* 9 */ "basefreq", /* 10 */ "rates", /* 11 */ "dimensions", /* 12 */ "ntax", /* 13 */ "nchar", /* 14 */ "format", /* 15 */ "datatype", /* 16 */ "yes", /* 17 */ "no", /* 18 */ "dna", /* 19 */ "rna", /* 20 */ "protein", /* 21 */ "interleave", /* 22 */ "gap", /* 23 */ "matchchar", /* 24 */ "missing", /* 25 */ "matrix", /* 26 */ "nexus", /* 27 */ "equal", /* 28 */ "gamma", /* 29 */ "sitespec", /* 30 */ "covarion", /* 31 */ "clock", /* 32 */ "unconstrained", /* 33 */ "relaxed", /* 34 */ "strict", /* 35 */ "empirical", /* 36 */ "estimate", /* 37 */ "tratio", /* 38 */ "revmat", /* 39 */ "nonrevmat", /* 40 */ "shape", /* 41 */ "sitepartition", /* 42 */ "siterates", /* 43 */ "partition", /* 44 */ "charset", /* 45 */ "exclude", /* 46 */ "outgroup", /* 47 */ "charstat", /* 48 */ "calibration", /* 49 */ "constraint", /* 50 */ "compchar", /* 51 */ "all", /* 52 */ "include", /* 53 */ "delete", /* 54 */ "restore", /* 55 */ "mcmc", /* 56 */ "seed", /* 57 */ "ngen", /* 58 */ "samplefreq", /* 59 */ "printfreq", /* 60 */ "nchains", /* 61 */ "temp", /* 62 */ "filename", /* 63 */ "prset", /* 64 */ "qmatpr", /* 65 */ "fixed", /* 66 */ "exp", /* 67 */ "uni", /* 68 */ "basefreqpr", /* 69 */ "dir", /* 70 */ "brlenpr", /* 71 */ "bd", /* 72 */ "shapepr", /* 73 */ "siteratepr", /* 74 */ "startingtree", /* 75 */ "random", /* 76 */ "user", /* 77 */ "usertree", /* 78 */ "showtree", /* 79 */ "enforcecal", /* 80 */ "enforcecon", /* 81 */ "inferanc", /* 82 */ "root", /* 83 */ "deroot", /* 84 */ "ncat", /* 85 */ "shownodes", /* 86 */ "savebrlens", /* 87 */ "seqerror", /* 88 */ "sumt", /* 89 */ "sump", /* 90 */ "burnin", /* 91 */ "tree", /* 92 */ "translate", /* 93 */ "contype", /* 94 */ "halfcompat", /* 95 */ "allcompat", /* 96 */ "bayesparams", /* 97 */ "header", /* 98 */ "set", /* 99 */ "autoclose", /* 100 */ "mcmcp", /* 101 */ "nperbs", /* 102 */ "aamodel", /* 103 */ "poisson", /* 104 */ "equalin", /* 105 */ "jones", /* 106 */ "dayhoff", /* 107 */ "gtr", /* 108 */ "omegapr", /* 109 */ "ancfile", /* 110 */ "enforcecodon", /* 111 */ "omegavar", /* 112 */ "ny98", /* 113 */ "omega", /* 114 */ "code", /* 115 */ "universal", /* 116 */ "vertmt", /* 117 */ "mycoplasma", /* 118 */ "yeast", /* 119 */ "ciliates", /* 120 */ "metmt", /* 121 */ "inferpossel", /* 122 */ "adgamma", /* 123 */ "ac1ny98", /* 124 */ "ac2ny98", /* 125 */ "lag", /* 126 */ "lagpr", /* 127 */ "invariant", /* 128 */ "invgamma", /* 129 */ "ssgamma", /* 130 */ "ssadgamma", /* 131 */ "inferrates", /* 132 */ "covny98", /* 133 */ "restriction", /* 134 */ "qmatwin", /* 135 */ "pialpha", /* 136 */ "ratewin", /* 137 */ "tuning", /* 138 */ "clocktune", /* 139 */ "timetune", /* 140 */ "mtune", /* 141 */ "extendp", /* 142 */ "wormtune", /* 143 */ "bdwin", /* 144 */ "omegawin", /* 145 */ "omegaalpha", /* 146 */ "rhowin", /* 147 */ "invslide", /* 148 */ "switchwin", /* 149 */ "alphawin", /* 150 */ "props", /* 151 */ "reconlim", /* 152 */ "covtype", /* 153 */ "ts98", /* 154 */ "gcswitch", /* 155 */ "tbrextendp", /* 156 */ "tbrtune", /* 157 */ "bltune", /* 158 */ "standard", /* 159 */ "mixed", /* 160 */ "coding", /* 161 */ "variable", /* 162 */ "informative", /* 163 */ "noabsence", /* 164 */ "nopresence", /* 165 */ "statefreqpr", /* 166 */ "setmoveprops", /* 167 */ "cyclefreq", /* 168 */ "calcpseudobf", /* 169 */ "parsmodel", /* 170 */ "typeset", /* 171 */ "nstatesset", /* 172 */ "ord", /* 173 */ "unord", /* 174 */ "assume"}; /* 175 */ #if defined(__MWERKS__) # pragma optimization_level 1 #endif char ChangeCase (char c) { int x; x = tolower(c); return (x); } void CleanExit (void) { if (isMatrixAllocated == YES) { free (origDataMatrix); free (taxonLabels); free (excludedSites); free (sitePartitions); free (matrixCounter); free (charSet); free (calibrations); free (constraints); free (excludedTaxa); free (userTree); free (sumTree); free (transOldNames); free (transNewNames); free (charType); free (charNStates); free (typeDef); free (nStatesDef); } nTaxa = nChar = 0; isMatrixAllocated = NO; isNtaxaInitialized = NO; isNcharInitialized = NO; } int Dex(command_list *p) { if (p == NULL) return (-1); else return (p->index); } int Dex2 (TreeNode *p) { return (p == NULL) ? -1 : p->index; } int Dex3 (ConTreeNode *p) { return (p == NULL) ? -1 : p->index; } int DerootTree (TreeNode *p) { int i, nNodes, isMarked; double tempBrLen; TreeNode **downPass, *lft, *m1, *m2, *um1, *um2, *out; if (isTreeRooted == NO) { printf (" ERROR: Tree is already unrooted\n"); return (ERROR); } else nNodes = 2 * nTaxa; /* bring the outgroup around to the first right position */ downPass = (TreeNode **)malloc((size_t) (2 * nTaxa * sizeof(TreeNode *))); if (!downPass) { printf (" ERROR: Could not allocate downPass\n"); return (ERROR); } i = 0; GetUserDownPass (userTreeRoot, downPass, &i); isMarked = NO; for (i=0; iindex == outgroupNum) { p->marked = YES; isMarked = YES; } else p->marked = NO; } if (isMarked == NO) { printf (" ERROR: Could not find outgroup taxon\n"); return (ERROR); } for (i=0; ileft != NULL && p->right != NULL) if (p->left->marked == YES || p->right->marked == YES) p->marked = YES; } free (downPass); lft = userTreeRoot->left; while (lft->left->index != outgroupNum && lft->right->index != outgroupNum) { 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 = lft; m1->left = um2; m1->right = um1; m1->anc = lft; um1->anc = m1; um2->anc = m1; m1->marked = NO; } else { printf (" ERROR: Rooting routine is lost\n"); return (ERROR); } } 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 = lft; m1->left = um1; m1->right = um2; m1->anc = lft; um1->anc = m1; um2->anc = m1; m1->marked = NO; } else { printf (" ERROR: Rooting routine is lost\n"); return (ERROR); } } else { printf (" ERROR: Rooting routine is lost\n"); return (ERROR); } } /* make certain outgroup is to the right of the root */ if (userTreeRoot->left->left->index == outgroupNum) { m1 = userTreeRoot->left->left; m2 = userTreeRoot->left->right; lft = userTreeRoot->left; lft->left = m2; lft->right = m1; } /* now, take outgroup and make it point down */ m1 = userTreeRoot; m2 = userTreeRoot->left; lft = userTreeRoot->left->left; out = userTreeRoot->left->right; tempBrLen = out->length; if (tempBrLen < 0.0) tempBrLen = 0.0; lft->anc = out; out->left = lft; out->right = NULL; out->anc = NULL; lft->length += tempBrLen; m1->left = NULL; m2->left = NULL; m1->right = NULL; m2->right = NULL; m1->anc = NULL; m2->anc = NULL; userTreeRoot = out; isTreeRooted = NO; /* reindex internal nodes of tree */ i = nTaxa; FinishTree (userTreeRoot, &i, NO); return (NO_ERROR); } int DerootSumTree (TreeNode *p) { int i, nNodes, isMarked; double tempBrLen; TreeNode **downPass, *lft, *m1, *m2, *um1, *um2, *out; if (isSumTreeRooted == NO) { printf (" ERROR: Tree is already unrooted\n"); return (ERROR); } else nNodes = 2 * nTaxa; /* bring the outgroup around to the first right position */ downPass = (TreeNode **)malloc((size_t) (2 * nTaxa * sizeof(TreeNode *))); if (!downPass) { printf (" ERROR: Could not allocate downPass\n"); return (ERROR); } i = 0; GetUserDownPass (sumTreeRoot, downPass, &i); isMarked = NO; for (i=0; iindex == outgroupNum) { p->marked = YES; isMarked = YES; } else p->marked = NO; } if (isMarked == NO) { printf (" ERROR: Could not find outgroup taxon\n"); return (ERROR); } for (i=0; ileft != NULL && p->right != NULL) if (p->left->marked == YES || p->right->marked == YES) p->marked = YES; } free (downPass); lft = sumTreeRoot->left; while (lft->left->index != outgroupNum && lft->right->index != outgroupNum) { 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 = lft; m1->left = um2; m1->right = um1; m1->anc = lft; um1->anc = m1; um2->anc = m1; m1->marked = NO; } else { printf (" ERROR: Rooting routine is lost\n"); return (ERROR); } } 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 = lft; m1->left = um1; m1->right = um2; m1->anc = lft; um1->anc = m1; um2->anc = m1; m1->marked = NO; } else { printf (" ERROR: Rooting routine is lost\n"); return (ERROR); } } else { printf (" ERROR: Rooting routine is lost\n"); return (ERROR); } } /* make certain outgroup is to the right of the root */ if (sumTreeRoot->left->left->index == outgroupNum) { m1 = sumTreeRoot->left->left; m2 = sumTreeRoot->left->right; lft = sumTreeRoot->left; lft->left = m2; lft->right = m1; } /* now, take outgroup and make it point down */ m1 = sumTreeRoot; m2 = sumTreeRoot->left; lft = sumTreeRoot->left->left; out = sumTreeRoot->left->right; tempBrLen = out->length; if (tempBrLen < 0.0) tempBrLen = 0.0; lft->anc = out; out->left = lft; out->right = NULL; out->anc = NULL; lft->length += tempBrLen; m1->left = NULL; m2->left = NULL; m1->right = NULL; m2->right = NULL; m1->anc = NULL; m2->anc = NULL; sumTreeRoot = out; isSumTreeRooted = NO; /* reindex internal nodes of tree */ i = nTaxa; FinishTree (sumTreeRoot, &i, NO); return (NO_ERROR); } void DisplayBits (unsigned value) { int unsignedIntSize; unsigned c, displayMask; unsignedIntSize = sizeof (unsigned); displayMask = 1 << (8*unsignedIntSize-1); //printf ("%10u = ", value); for (c=1; c<=8*unsignedIntSize; c++) { putchar (value & displayMask ? '1' : '0'); value <<= 1; if (c % 8 == 0) putchar (' '); } //putchar ('\n'); } int DoAction (int actionType) { int i, j, k, l, nuc, tempInt, tempInt2, foundSet, nNodes, offSet, foundPartition, allSame, nextConNode, isCompat, whichIsX, whichIsY, numY[SCREENWIDTH], unsignedIntsize; unsigned x, y, lastMask; long int tempLongInt; double tempD, *tempInfo, a, aOld, s, minX, minY, maxX, maxY, meanY[SCREENWIDTH]; char tempString[100], tempString2[100], conFile[100]; TreeNode *p; FILE *fp; # if defined (CONTREE) PolyNode *conRoot, *cp, *q, *r, *ql, *pl, *conNodes; int targetNode, nBits; double freq; unsigned int *partition, *ingroupPartition, *outgroupPartition; # else ConTreeNode *conRoot, *cp, *cl, *cr, *cn, *cq, **taxaInConstraint, **taxaNotInConstraint; unsigned *compatParts, *tempParts; int nTaxaInConstraint, nTaxaNotInConstraint, nCompatParts, nSame1, nSame2, *numOnes; double *partLengths; # endif switch (actionType) { case CMD_QUIT: printf (" Quitting program\n\n"); CleanExit(); exit (0); break; case CMD_EXECUTE: printf (" Executing file named \"%s\"\n", inFile); if (ReadNexusFile () == ERROR) { printf (" Could not execute file\n"); ErrorExit (); return (ERROR); } break; case CMD_ENTER_INPUT_FILE_NAME1: changeToLower = NO; break; case CMD_ENTER_INPUT_FILE_NAME2: changeToLower = YES; strcpy (inFile, token); strcpy (outFile, inFile); strcat (outFile, ".out"); strcpy (sumTreeFile, token); strcat (sumTreeFile, ".out.t"); strcpy (sumParmFile, token); strcat (sumParmFile, ".out.bp"); break; case CMD_INIT_NTAX: sscanf (token, "%d", &tempInt); nTaxa = tempInt; isNtaxaInitialized = YES; break; case CMD_INIT_NCHAR: sscanf (token, "%d", &tempInt); nChar = tempInt; isNcharInitialized = YES; break; case CMD_INIT_DATA_BLOCK: if (isMatrixAllocated == YES) { free (origDataMatrix); free (taxonLabels); free (excludedSites); free (sitePartitions); free (matrixCounter); free (charSet); free (calibrations); free (constraints); free (excludedTaxa); free (userTree); free (charType); free (charNStates); free (typeDef); free (nStatesDef); nTaxa = nChar = 0; } isMatrixAllocated = NO; isNtaxaInitialized = NO; isNcharInitialized = NO; isUserTreeDefined = NO; isInterleaved = NO; gapChar = 0; matchChar = 0; missingChar = 0; dataType = -1; printf (" Entering data block\n"); break; case CMD_ALLOC_MATRIX: if (isMatrixAllocated == NO && isNtaxaInitialized == YES && isNcharInitialized == YES) { nCharStates = 4; origDataMatrix = (int *)malloc((size_t) (nTaxa * nChar * sizeof(int))); if (!origDataMatrix) { printf (" ERROR: Could not allocate origDataMatrix\n"); ErrorExit (); return (ERROR); } for (i=0; i= nTaxa) { if (isInterleaved == YES) whichTaxon = 0; else { printf (" ERROR: Too many taxa (line %d)\n", lineNumber); ErrorExit (); return (ERROR); } } if (isFirstInterleavedBlock == YES) { i = 0; while (token[i] != '\0') { taxonLabels[whichTaxon*100+i] = token[i]; i++; } taxonLabels[whichTaxon*100+i] = '|'; returnEncountered = NO; whichTaxon++; taxonLabels[nTaxa*100-1] = '\0'; } else { i = 0; while (taxonLabels[whichTaxon*100+i] != '|') { tempString[i] = taxonLabels[whichTaxon*100+i]; i++; } tempString[i] = '\0'; if (!strcmp(tempString, token)) { printf (" ERROR: Taxon names don't match across interleaved blocks (line %d)\n", lineNumber); ErrorExit (); return (ERROR); } } } else { i = 0; while (token[i] != '\0') { if (matrixCounter[whichTaxon-1] == nChar) { printf (" ERROR: Too many characters (line %d)\n", lineNumber); ErrorExit (); return (ERROR); } //printf ("%2d %2d -- %c\n", whichTaxon, matrixCounter[whichTaxon-1], token[i]); if (token[i] == matchChar) { if (whichTaxon - 1 == 0) { printf (" ERROR: Matching characters (%c) cannot be in first taxon (line %d)\n", matchChar, lineNumber); ErrorExit (); return (ERROR); } origDataMatrix[(whichTaxon-1)*nChar+matrixCounter[whichTaxon-1]] = origDataMatrix[(0)*nChar+matrixCounter[whichTaxon-1]]; matrixCounter[whichTaxon-1]++; } else { if (dataType == DNA || dataType == RNA) nuc = NucID (token[i]); else if (dataType == PROTEIN) nuc = ProtID (token[i]); else if (dataType == RESTRICTION) nuc = ResID (token[i]); else if (dataType == STANDARD) nuc = StandID (token[i]); else if (dataType == MIXED) { nuc = 0; printf (" ERROR: Unrecognized data type\n"); ErrorExit (); return (ERROR); } else { printf (" ERROR: Unrecognized data type\n"); ErrorExit (); return (ERROR); } if (nuc != -1) { origDataMatrix[(whichTaxon-1)*nChar+matrixCounter[whichTaxon-1]] = nuc; matrixCounter[whichTaxon-1]++; } else { printf (" ERROR: Unrecognized character state %c (line %d)\n", token[i], lineNumber); ErrorExit (); return (ERROR); } } i++; } } } else { printf (" ERROR: Cannot read in matrix because it is not allocated\n"); } break; case CMD_FINISH_READ_MATRIX: for (i=0; i "); for (i=0; i 5) { printf (" ERROR: Too many user rates assigned\n"); return (ERROR); } else { sscanf (token, "%lf", &tempD); revRates[globalInt++] = tempD; } break; case CMD_INIT_REVRATES3: if (globalInt != 6) { printf (" ERROR: Incorrect number of rates assigned\n"); return (ERROR); } else { for (i=0; i<5; i++) revRates[i] /= revRates[5]; revRates[5] = 1.0; } break; case CMD_INIT_NONREVRATES: if (token[0] == 'e') { strcpy(nonRevRatesModel, token); tRatio = 1.0; } else { globalInt = 0; strcpy (nonRevRatesModel, "user defined"); } printf (" Setting rates of non-reversible model\n", tRatio); break; case CMD_INIT_NONREVRATES2: if (globalInt > 11) { printf (" ERROR: Too many user rates assigned\n"); return (ERROR); } else { sscanf (token, "%lf", &tempD); nonRevRates[globalInt++] = tempD; } break; case CMD_INIT_NONREVRATES3: if (globalInt != 12) { printf (" ERROR: Incorrect number of rates assigned\n"); return (ERROR); } else { for (i=0; i<11; i++) nonRevRates[i] /= nonRevRates[11]; nonRevRates[11] = 1.0; } break; case CMD_INIT_SHAPE: if (token[0] == 'e') { strcpy(shapeModel, token); shape = 0.5; } else { sscanf (token, "%lf", &tempD); shape = tempD; strcpy (shapeModel, "user defined"); } printf (" Setting gamma shape parameter = %lf\n", shape); break; case CMD_INIT_PARTITION: strcpy(partitionName, token); break; case CMD_INIT_SITERATES: if (token[0] == 'e') { strcpy(siteRatesModel, token); } else { globalInt = 0; strcpy (siteRatesModel, "user defined"); } printf (" Setting rates for partitions\n"); break; case CMD_INIT_SITERATES2: if (globalInt > 49) { printf (" ERROR: Too many site rates assigned\n"); return (ERROR); } else { sscanf (token, "%lf", &tempD); siteRates[globalInt++] = tempD; } break; case CMD_INIT_CHARSET_NAME: if (isMatrixAllocated == NO) { printf (" ERROR: Data not read in\n"); return (ERROR); } i = 0; while (token[i] != '\0') { charSetLabels[nCharSets][i] = token[i]; i++; } charSetLabels[nCharSets][i] = '|'; charSetFrom = -1; charSetTo = -1; charSetWait = NO; charSetSlashWait = NO; charSetPotWait = YES; charSetFreq = 1; printf (" Assigning character set \""); i = 0; while (charSetLabels[nCharSets][i] != '|') { printf ("%c", charSetLabels[nCharSets][i]); i++; } printf ("\"\n"); nCharSets++; break; case CMD_INIT_CHARSET_NAME2: if (isMatrixAllocated == NO) { printf (" ERROR: Data not read in\n"); return (ERROR); } printf (" ERROR: You cannot use charset names that are digits\n"); return (ERROR); break; case CMD_INIT_CHARSET: if (token[0] == '-') { charSetWait = YES; charSetPotWait = YES; } else if (token[0] == '\\') { charSetSlashWait = YES; charSetPotWait = YES; } else { if (charSetWait == NO && charSetSlashWait == NO && charSetPotWait == NO) { j = charSetFreq; for (i=charSetFrom-1; i<=charSetTo-1; i++) { if (j % charSetFreq == 0) charSet[i] += (int)pow(2.0, (double)nCharSets); j++; } charSetFrom = -1; charSetTo = -1; charSetWait = NO; charSetSlashWait = NO; charSetFreq = 1; charSetPotWait = YES; } sscanf (token, "%d", &tempInt); if (charSetSlashWait == YES && charSetWait == YES) { printf (" ERROR: Problem assigning character set\n"); return (ERROR); } if (charSetSlashWait == YES) { charSetFreq = tempInt; charSetSlashWait = NO; charSetPotWait = NO; } else if (charSetWait == YES) { charSetTo = tempInt; charSetWait = NO; charSetPotWait = NO; } else { if (charSetWait == NO) { charSetFrom = tempInt; charSetTo = tempInt; charSetPotWait = NO; } else { charSetTo = tempInt; charSetPotWait = NO; } } } break; case CMD_INIT_CHARSET2: if (charSetWait == NO && charSetSlashWait == NO && charSetPotWait == NO) { j = charSetFreq; for (i=charSetFrom-1; i<=charSetTo-1; i++) { if (j % charSetFreq == 0) charSet[i] += (int)pow(2.0, (double)nCharSets); j++; } charSetFrom = -1; charSetTo = -1; charSetWait = NO; charSetSlashWait = NO; charSetFreq = 1; charSetPotWait = YES; } # if 0 for (i=0; i nPartitions) { printf (" ERROR: You exceeded the %d partitions you specified.\n", nPartitions); return (ERROR); } break; case CMD_PARTITION4: if (token[0] == '-') { charSetWait = YES; charSetPotWait = YES; } else if (token[0] == '\\') { charSetSlashWait = YES; charSetPotWait = YES; } else { if (charSetWait == NO && charSetSlashWait == NO && charSetPotWait == NO) { j = charSetFreq; for (i=charSetFrom-1; i<=charSetTo-1; i++) { if (j % charSetFreq == 0) sitePartitions[i] = whichPartition; j++; } charSetFrom = -1; charSetTo = -1; charSetWait = NO; charSetSlashWait = NO; charSetFreq = 1; charSetPotWait = YES; } sscanf (token, "%d", &tempInt); if (charSetSlashWait == YES && charSetWait == YES) { printf (" ERROR: Problem assigning character set\n"); return (ERROR); } if (charSetSlashWait == YES) { charSetFreq = tempInt; charSetSlashWait = NO; charSetPotWait = NO; } else if (charSetWait == YES) { charSetTo = tempInt; charSetWait = NO; charSetPotWait = NO; } else { if (charSetWait == NO) { charSetFrom = tempInt; charSetTo = tempInt; charSetPotWait = NO; } else { charSetTo = tempInt; charSetPotWait = NO; } } } break; case CMD_PARTITION5: if (charSetWait == NO && charSetSlashWait == NO && charSetPotWait == NO) { j = charSetFreq; for (i=charSetFrom-1; i<=charSetTo-1; i++) { if (j % charSetFreq == 0) sitePartitions[i] = whichPartition; j++; } charSetFrom = -1; charSetTo = -1; charSetWait = NO; charSetSlashWait = NO; charSetFreq = 1; charSetPotWait = YES; } if (token[0] == ';') { if (whichPartition != nPartitions) { printf (" ERROR: An incorrect number of partitions are defined.\n", nPartitions); return (ERROR); } isPartitionDefined = YES; } if (token[0] == ',') { whichPartition++; if (whichPartition > nPartitions) { printf (" ERROR: You exceeded the %d partitions you specified.\n", nPartitions); return (ERROR); } } else { for (i=0; i nTaxa) { printf (" ERROR: Outgroup number out of range\n"); return (ERROR); } else { tempInt--; outgroupNum = tempInt; printf (" Setting outgroup to species \""); i = 0; while (taxonLabels[outgroupNum*100+i] != '|') { printf ("%c", taxonLabels[outgroupNum*100+i]); i++; } printf ("\"\n"); } if (isUserTreeDefined == YES) { if (isTreeRooted == YES) { if (DerootTree (userTreeRoot) == ERROR) return (ERROR); if (RootTree (userTreeRoot) == ERROR) return (ERROR); } else { if (RootTree (userTreeRoot) == ERROR) return (ERROR); if (DerootTree (userTreeRoot) == ERROR) return (ERROR); } } break; case CMD_SHOW_CHAR_STAT: if (isMatrixAllocated == NO) { printf (" ERROR: Data not read in\n"); return (ERROR); } printf (" site included/excluded partition number\n"); printf (" ------------------------------------------------\n"); for (i=0; i= 31) { printf (" ERROR: Too many calibrations (32 maximum)\n"); return (ERROR); } printf (" Reading calibration \"%s\"\n", token); i = 0; while (token[i] != '\0') { calibrationLabels[nCalibrations][i] = token[i]; i++; } calibrationLabels[nCalibrations][i] = '|'; break; case CMD_CALIBRATION2: printf (" ERROR: Calibration label cannot be a number\n"); break; case CMD_CALIBRATION3: if (isMatrixAllocated == NO) { printf (" ERROR: Data not read in\n"); return (ERROR); } if (nCalibrations >= 31) { printf (" ERROR: Too many calibrations (32 maximum)\n"); return (ERROR); } sscanf (token, "%lf", &tempD); calibrationTimes[nCalibrations][0] = tempD; break; case CMD_CALIBRATION4: if (isMatrixAllocated == NO) { printf (" ERROR: Data not read in\n"); return (ERROR); } if (nCalibrations >= 31) { printf (" ERROR: Too many calibrations (32 maximum)\n"); return (ERROR); } sscanf (token, "%lf", &tempD); calibrationTimes[nCalibrations][1] = tempD; break; case CMD_CALIBRATION5: if (isMatrixAllocated == NO) { printf (" ERROR: Data not read in\n"); return (ERROR); } if (nCalibrations >= 31) { printf (" ERROR: Too many calibrations (32 maximum)\n"); return (ERROR); } tempInt = NO; for (j=0; j= 31) { printf (" ERROR: Too many calibrations (32 maximum)\n"); return (ERROR); } sscanf (token, "%d", &tempInt); tempInt--; if (tempInt < 0 || tempInt >= nTaxa) { printf (" ERROR: Taxon number is out of range\n"); } else { calibrations[tempInt] += (int)pow(2.0, (double)(nCalibrations)); } break; case CMD_CALIBRATION7: # if 0 i = 0; while (calibrationLabels[nCalibrations][i] != '|') { printf ("%c", calibrationLabels[nCalibrations][i]); i++; } printf (" (%lf %lf) -> ", calibrationTimes[nCalibrations][0], calibrationTimes[nCalibrations][1]); for (i=0; i= 31) { printf (" ERROR: Too many constraints (32 maximum)\n"); return (ERROR); } printf (" Reading constraint \"%s\"\n", token); i = 0; while (token[i] != '\0') { constraintLabels[nConstraints][i] = token[i]; i++; } constraintLabels[nConstraints][i] = '|'; break; case CMD_CONSTRAINT2: if (isMatrixAllocated == NO) { printf (" ERROR: Data not read in\n"); return (ERROR); } if (nConstraints >= 31) { printf (" ERROR: Too many constraints (32 maximum)\n"); return (ERROR); } printf (" ERROR: Constraint label cannot be a number\n"); break; case CMD_CONSTRAINT3: if (isMatrixAllocated == NO) { printf (" ERROR: Data not read in\n"); return (ERROR); } if (nConstraints >= 31) { printf (" ERROR: Too many constraints (32 maximum)\n"); return (ERROR); } tempInt = NO; for (j=0; j= 31) { printf (" ERROR: Too many constraints (32 maximum)\n"); return (ERROR); } sscanf (token, "%d", &tempInt); tempInt--; if (tempInt < 0 || tempInt >= nTaxa) { printf (" ERROR: Taxon number is out of range\n"); } else { constraints[tempInt] += (int)pow(2.0, (double)(nConstraints)); } break; case CMD_CONSTRAINT5: if (isMatrixAllocated == NO) { printf (" ERROR: Data not read in\n"); return (ERROR); } if (nConstraints >= 31) { printf (" ERROR: Too many constraints (32 maximum)\n"); return (ERROR); } # if 0 i = 0; while (constraintLabels[nConstraints][i] != '|') { printf ("%c", constraintLabels[nConstraints][i]); i++; } printf (" -> "); for (i=0; i= nTaxa) { printf (" ERROR: Taxon number is out of range\n"); } else { excludedTaxa[tempInt] = YES; printf (" Excluded taxon "); i = 0; while (taxonLabels[tempInt*100+i] != '|') { printf ("%c", taxonLabels[tempInt*100+i]); i++; } printf ("\n"); } break; case CMD_DELETE3: if (isMatrixAllocated == NO) { printf (" ERROR: Data not read in\n"); return (ERROR); } printf (" Excluded all taxa\n", token); for (i=0; i= nTaxa) { printf (" ERROR: Taxon number is out of range\n"); } else { excludedTaxa[tempInt] = NO; printf (" Included taxon "); i = 0; while (taxonLabels[tempInt*100+i] != '|') { printf ("%c", taxonLabels[tempInt*100+i]); i++; } printf ("\n"); } break; case CMD_RESTORE3: if (isMatrixAllocated == NO) { printf (" ERROR: Data not read in\n"); return (ERROR); } printf (" Included all taxa\n", token); for (i=0; i= 2*nTaxa) { printf (" ERROR: Too many nodes in tree\n"); return (ERROR); } if (pPtr->left == NULL) { pPtr = &userTree[nextAvailableNode]; nextAvailableNode++; qPtr->left = pPtr; pPtr->anc = qPtr; qPtr = pPtr; } else if (pPtr->right == NULL) { pPtr = &userTree[nextAvailableNode]; nextAvailableNode++; qPtr->right = pPtr; pPtr->anc = qPtr; qPtr = pPtr; } else if (pPtr->anc == NULL) { pPtr = &userTree[nextAvailableNode]; nextAvailableNode++; qPtr->anc = pPtr; pPtr->left = qPtr; qPtr = pPtr; userTreeRoot = pPtr; pPtr->marked = YES; } else { printf ("\n ERROR: Tree is not bifurcating\n"); return (ERROR); } } break; case CMD_TREE2: /* */ if (nextAvailableNode+1 >= 2*nTaxa) { printf ("\n ERROR: Too many nodes in tree\n"); return (ERROR); } tempInt = NO; for (j=0; jleft == NULL) { pPtr = &userTree[nextAvailableNode]; strcpy (pPtr->label, token); pPtr->index = j; nextAvailableNode++; qPtr->left = pPtr; pPtr->anc = qPtr; qPtr = pPtr; } else if (pPtr->right == NULL) { pPtr = &userTree[nextAvailableNode]; strcpy (pPtr->label, token); pPtr->index = j; nextAvailableNode++; qPtr->right = pPtr; pPtr->anc = qPtr; qPtr = pPtr; } else if (pPtr->anc == NULL) { pPtr = &userTree[nextAvailableNode]; strcpy (pPtr->label, token); pPtr->index = j; nextAvailableNode++; qPtr->anc = pPtr; pPtr->left = qPtr; qPtr = pPtr; userTreeRoot = pPtr; pPtr->marked = YES; } else { printf ("\n ERROR: Tree is not bifurcating\n"); return (ERROR); } } else { printf ("\n ERROR: Taxon name \"%s\" is not found\n", token); return (ERROR); } break; case CMD_TREE3: /* */ if (nextAvailableNode+1 >= 2*nTaxa) { printf ("\n ERROR: Too many nodes in tree\n"); return (ERROR); } sscanf (token, "%d", &tempInt); tempInt--; if (tempInt < 0 || tempInt >= nTaxa) { printf ("\n ERROR: Taxon number is out of range\n"); } else { if (pPtr->left == NULL) { pPtr = &userTree[nextAvailableNode]; i = 0; while (taxonLabels[tempInt*100+i] != '|') { pPtr->label[i] = taxonLabels[tempInt*100+i]; i++; } pPtr->label[i] = '\0'; pPtr->index = tempInt; nextAvailableNode++; qPtr->left = pPtr; pPtr->anc = qPtr; qPtr = pPtr; } else if (pPtr->right == NULL) { pPtr = &userTree[nextAvailableNode]; i = 0; while (taxonLabels[tempInt*100+i] != '|') { pPtr->label[i] = taxonLabels[tempInt*100+i]; i++; } pPtr->label[i] = '\0'; pPtr->index = tempInt; nextAvailableNode++; qPtr->right = pPtr; pPtr->anc = qPtr; qPtr = pPtr; } else if (pPtr->anc == NULL) { pPtr = &userTree[nextAvailableNode]; i = 0; while (taxonLabels[tempInt*100+i] != '|') { pPtr->label[i] = taxonLabels[tempInt*100+i]; i++; } pPtr->label[i] = '\0'; pPtr->index = tempInt; nextAvailableNode++; qPtr->anc = pPtr; pPtr->left = qPtr; qPtr = pPtr; userTreeRoot = pPtr; pPtr->marked = YES; } else { printf ("\n ERROR: Tree is not bifurcating\n"); return (ERROR); } } break; case CMD_TREE4: /* ) */ if (pPtr->marked == NO) { if (pPtr->anc != NULL) { pPtr = pPtr->anc; qPtr = pPtr; } else { printf ("\n ERROR: Cannot go down\n"); return (ERROR); } } else { if (pPtr->left != NULL) { pPtr = pPtr->left; qPtr = pPtr; } else { printf ("\n ERROR: Cannot go up\n"); return (ERROR); } } break; case CMD_TREE5: /* , */ if (pPtr->marked == NO) { if (pPtr->anc != NULL) { pPtr = pPtr->anc; qPtr = pPtr; } else { printf ("\n ERROR: Cannot go down\n"); return (ERROR); } } else { if (pPtr->left != NULL) { pPtr = pPtr->left; qPtr = pPtr; } else { printf ("\n ERROR: Cannot go up\n"); return (ERROR); } } break; case CMD_TREE6: /* : */ break; case CMD_TREE7: /* ; */ if (pPtr != &userTree[0]) { printf ("\n ERROR: Something wrong with user tree\n"); return (ERROR); } if (pPtr->anc == NULL) isTreeRooted = YES; else isTreeRooted = NO; if (isTreeRooted == YES) { printf (" User tree defined (rooted)\n"); pPtr = &userTree[nextAvailableNode]; nextAvailableNode++; pPtr->left = qPtr; qPtr->anc = pPtr; userTreeRoot = pPtr; } else { printf (" User tree defined (unrooted)\n"); } i = nTaxa; if (isTreeRooted == YES) FinishTree (userTreeRoot, &i, YES); else FinishTree (userTreeRoot, &i, NO); isUserTreeDefined = YES; if (isTreeRooted == YES) { if (DerootTree (userTreeRoot) == ERROR) return (ERROR); if (RootTree (userTreeRoot) == ERROR) return (ERROR); } else { if (RootTree (userTreeRoot) == ERROR) return (ERROR); if (DerootTree (userTreeRoot) == ERROR) return (ERROR); } break; case CMD_TREE8: /* */ sscanf (token, "%lf", &tempD); if (pPtr->marked == NO) pPtr->length = tempD; else { if (pPtr->left != NULL) pPtr->left->length = tempD; else { printf (" ERROR: Cannot assign branch length to left node\n"); printf ("token = %s\n", token); return (ERROR); } } userBrlensDef = YES; break; case CMD_SHOWUSERTREE: if (isUserTreeDefined == NO) { printf (" ERROR: No user tree defined\n"); } else { if (isTreeRooted == YES) printf ("\n Tree is rooted:\n\n"); else printf ("\n Tree is unrooted:\n\n"); if (ShowTree (userTreeRoot, isTreeRooted) == ERROR) { return (ERROR); } else printf ("\n"); } break; case CMD_INIT_ENFORCECAL: if (token[0] == 'y') enforceCalibrations = YES; else enforceCalibrations = NO; printf (" Setting enforcecal=%s\n", token); break; case CMD_INIT_ENFORCECON: if (token[0] == 'y') enforceConstraints = YES; else enforceConstraints = NO; printf (" Setting enforcecon=%s\n", token); break; case CMD_INIT_INFERANC: if (token[0] == 'y') inferAncStates = YES; else inferAncStates = NO; printf (" Setting inferanc=%s\n", token); break; case CMD_INIT_NCAT: sscanf (token, "%d", &tempInt); nRateCats = tempInt; printf (" Setting ncat=%d\n", nRateCats); break; case CMD_ROOT_TREE: if (isUserTreeDefined == NO) { printf (" ERROR: No user tree defined\n"); } else { if (isTreeRooted == YES) { printf ("\n WARNING: Tree is already rooted\n"); if (DerootTree (userTreeRoot) == ERROR) return (ERROR); if (RootTree (userTreeRoot) == ERROR) return (ERROR); } else { printf (" Rooting tree\n"); if (RootTree (userTreeRoot) == ERROR) return (ERROR); } } break; case CMD_DEROOT_TREE: if (isUserTreeDefined == NO) { printf (" ERROR: No user tree defined\n"); } else { if (isTreeRooted == YES) { printf (" Derooting tree\n"); if (DerootTree (userTreeRoot) == ERROR) return (ERROR); } else { printf (" WARNING: Tree is already unrooted\n"); if (RootTree (userTreeRoot) == ERROR) return (ERROR); if (DerootTree (userTreeRoot) == ERROR) return (ERROR); } } break; case CMD_SHOWNODES: if (isUserTreeDefined == NO) { printf (" ERROR: No user tree defined\n"); } else { if (isTreeRooted == YES) printf ("\n Tree is rooted:\n\n"); else printf ("\n Tree is unrooted:\n\n"); ShowNodes (userTreeRoot, 2, isTreeRooted); } break; case CMD_SAVEBRLENS: if (token[0] == 'y') saveBrlens = YES; else saveBrlens = NO; break; case CMD_SET_SEQ_ERROR: sscanf (token, "%lf", &tempD); if (tempD < 0.0 || tempD > 1.0) { printf (" ERROR: Error probability must be between 0 and 1\n"); return (ERROR); } seqErrorProb = tempD; printf (" Setting sequencing error probability\n"); break; case CMD_SET_SP_RATE: sscanf (token, "%lf", &tempD); if (tempD < 0.0) { printf (" ERROR: Speciation rate must be greater than zero\n"); return (ERROR); } speciationRate = tempD; printf (" Setting speciation rate\n"); break; case CMD_SET_EX_RATE: sscanf (token, "%lf", &tempD); if (tempD < 0.0) { printf (" ERROR: Extinction rate must be greater than zero\n"); return (ERROR); } extinctionRate = tempD; printf (" Setting extinction rate\n"); break; case CMD_SET_SAMP_FRAC: sscanf (token, "%lf", &tempD); if (tempD < 0.0) { printf (" ERROR: Taxon sampling fraction must be between 0 and 1.\n"); return (ERROR); } samplingFraction = tempD; printf (" Setting taxon sampling fraction\n"); break; case CMD_SUMT_FILENAME: strcpy (sumTreeFile, token); break; case CMD_SUMT_BURN: sscanf (token, "%d", &tempInt); burnIn = tempInt; break; case CMD_SUMT_READFILE: if (isMatrixAllocated == NO) { printf (" ERROR: Data not read in\n"); return (ERROR); } else { foundBrlens = NO; if (SummarizeTrees () == ERROR) { printf (" Could not execute file\n"); ErrorExit (); return (ERROR); } } break; case CMD_TRANS_INIT: transNumber = 0; //printf (" Reading translation table\n"); break; case CMD_TRANS_NEW_NM: //printf (" %s", token); i = 0; while (token[i] != '\0') { transNewNames[transNumber*100+i] = token[i]; i++; } transNewNames[transNumber*100+i] = '|'; break; case CMD_TRANS_OLD_NM: //printf (" = %s\n", token); i = 0; while (token[i] != '\0') { transOldNames[transNumber*100+i] = token[i]; i++; } transOldNames[transNumber*100+i] = '|'; break; case CMD_TRANS_COMMA: transNumber++; break; case CMD_TRANS_END: transNumber++; /*for (j=0; j "); i = 0; while (transNewNames[j*100+i] != '|') { printf ("%c", transNewNames[j*100+i]); i++; } printf ("\n"); }*/ isTranslated = YES; break; case CMD_SUM_TREE: /* = */ if (isMatrixAllocated == NO) { printf (" ERROR: Data not read in\n"); return (ERROR); } for (i=0; i<2*nTaxa; i++) { sumTree[i].left = NULL; sumTree[i].right = NULL; sumTree[i].anc = NULL; sumTree[i].memoryIndex = i; sumTree[i].index = i; sumTree[i].upDateCl = NO; sumTree[i].length = 0.0; sumTree[i].marked = NO; strcpy (sumTree[i].label, "no name"); strcpy (sumTree[i].calibrationName, "no name"); strcpy (sumTree[i].constraintName, "no name"); } isFirstNode = YES; pPtr = qPtr = &sumTree[0]; nextAvailableNode = 0; break; case CMD_SUM_TREE1: /* ( */ if (isFirstNode == YES) { pPtr = &sumTree[nextAvailableNode]; nextAvailableNode++; isFirstNode = NO; sumTreeRoot = pPtr; } else { if (nextAvailableNode+1 >= 2*nTaxa) { printf (" ERROR: Too many nodes in tree\n"); return (ERROR); } if (pPtr->left == NULL) { pPtr = &sumTree[nextAvailableNode]; nextAvailableNode++; qPtr->left = pPtr; pPtr->anc = qPtr; qPtr = pPtr; } else if (pPtr->right == NULL) { pPtr = &sumTree[nextAvailableNode]; nextAvailableNode++; qPtr->right = pPtr; pPtr->anc = qPtr; qPtr = pPtr; } else if (pPtr->anc == NULL) { pPtr = &sumTree[nextAvailableNode]; nextAvailableNode++; qPtr->anc = pPtr; pPtr->left = qPtr; qPtr = pPtr; sumTreeRoot = pPtr; pPtr->marked = YES; } else { printf ("\n ERROR: Tree is not bifurcating\n"); return (ERROR); } } break; case CMD_SUM_TREE2: /* */ if (nextAvailableNode+1 >= 2*nTaxa) { printf ("\n ERROR: Too many nodes in tree\n"); return (ERROR); } if (isTranslated == YES) { /* need to use translate table */ /* first, find the name from the tree in the list of new labels */ tempInt = NO; for (j=0; jleft == NULL) { pPtr = &sumTree[nextAvailableNode]; strcpy (pPtr->label, tempString); pPtr->index = j; nextAvailableNode++; qPtr->left = pPtr; pPtr->anc = qPtr; qPtr = pPtr; } else if (pPtr->right == NULL) { pPtr = &sumTree[nextAvailableNode]; strcpy (pPtr->label, tempString); pPtr->index = j; nextAvailableNode++; qPtr->right = pPtr; pPtr->anc = qPtr; qPtr = pPtr; } else if (pPtr->anc == NULL) { pPtr = &sumTree[nextAvailableNode]; strcpy (pPtr->label, tempString); pPtr->index = j; nextAvailableNode++; qPtr->anc = pPtr; pPtr->left = qPtr; qPtr = pPtr; sumTreeRoot = pPtr; pPtr->marked = YES; } else { printf ("\n ERROR: Tree is not bifurcating\n"); return (ERROR); } break; case CMD_SUM_TREE3: /* */ if (nextAvailableNode+1 >= 2*nTaxa) { printf ("\n ERROR: Too many nodes in tree\n"); return (ERROR); } if (isTranslated == YES) { /* need to use translate table */ /* first, find the name from the tree in the list of new labels */ tempInt = NO; for (j=0; j= nTaxa) { printf ("\n ERROR: Taxon number is out of range\n"); } i = 0; while (taxonLabels[tempInt*100+i] != '|') { tempString[i] = taxonLabels[tempInt*100+i]; i++; } tempString[i] = '\0'; } /* add node */ if (pPtr->left == NULL) { pPtr = &sumTree[nextAvailableNode]; strcpy (pPtr->label, tempString); pPtr->index = tempInt; nextAvailableNode++; qPtr->left = pPtr; pPtr->anc = qPtr; qPtr = pPtr; } else if (pPtr->right == NULL) { pPtr = &sumTree[nextAvailableNode]; strcpy (pPtr->label, tempString); pPtr->index = tempInt; nextAvailableNode++; qPtr->right = pPtr; pPtr->anc = qPtr; qPtr = pPtr; } else if (pPtr->anc == NULL) { pPtr = &sumTree[nextAvailableNode]; strcpy (pPtr->label, tempString); pPtr->index = tempInt; nextAvailableNode++; qPtr->anc = pPtr; pPtr->left = qPtr; qPtr = pPtr; sumTreeRoot = pPtr; pPtr->marked = YES; } else { printf ("\n ERROR: Tree is not bifurcating\n"); return (ERROR); } break; case CMD_SUM_TREE4: /* ) */ if (pPtr->marked == NO) { if (pPtr->anc != NULL) { pPtr = pPtr->anc; qPtr = pPtr; } else { printf ("\n ERROR: Cannot go down\n"); return (ERROR); } } else { if (pPtr->left != NULL) { pPtr = pPtr->left; qPtr = pPtr; } else { printf ("\n ERROR: Cannot go up\n"); return (ERROR); } } break; case CMD_SUM_TREE5: /* , */ if (pPtr->marked == NO) { if (pPtr->anc != NULL) { pPtr = pPtr->anc; qPtr = pPtr; } else { printf ("\n ERROR: Cannot go down\n"); return (ERROR); } } else { if (pPtr->left != NULL) { pPtr = pPtr->left; qPtr = pPtr; } else { printf ("\n ERROR: Cannot go up\n"); return (ERROR); } } break; case CMD_SUM_TREE6: /* : */ break; case CMD_SUM_TREE7: /* ; */ numSumTrees++; if (numSumTrees % 1000 == 0) printf (" Reading tree %d\n", numSumTrees); //printf ("Reading tree %d\n", numSumTrees); if (pPtr != &sumTree[0]) { printf ("\n ERROR: Something wrong with user tree\n"); return (ERROR); } if (pPtr->anc == NULL) isSumTreeRooted = YES; else isSumTreeRooted = NO; if (isSumTreeRooted == YES) { //printf (" Tree read (rooted) %d\n", numSumTrees); pPtr = &sumTree[nextAvailableNode]; nextAvailableNode++; pPtr->left = qPtr; qPtr->anc = pPtr; sumTreeRoot = pPtr; } else { //printf (" Tree read (unrooted) %d\n", numSumTrees); } i = nTaxa; if (isSumTreeRooted == YES) FinishTree (sumTreeRoot, &i, YES); else FinishTree (sumTreeRoot, &i, NO); if (isSumTreeRooted == YES) { if (DerootSumTree (sumTreeRoot) == ERROR) return (ERROR); if (RootSumTree (sumTreeRoot) == ERROR) return (ERROR); } else { if (RootSumTree (sumTreeRoot) == ERROR) return (ERROR); if (DerootSumTree (sumTreeRoot) == ERROR) return (ERROR); } # if 0 if (ShowTree (sumTreeRoot, isSumTreeRooted) == ERROR) { printf (" ERROR: Could not show tree\n"); return (ERROR); } # endif /* allocate space for down pass sequence, if this is the first tree read in */ if (isSumTreeRooted == YES) nNodes = 2 * nTaxa; else nNodes = 2 * nTaxa - 2; if (numSumTrees == 1) { sumTreeDP = (TreeNode **)malloc((size_t) (2 * nTaxa * sizeof(TreeNode *))); if (!sumTreeDP) { printf (" ERROR: Could not allocate sumTreeDP\n"); return (ERROR); } } /* set up bits, if this is the first tree read in */ if (numSumTrees == 1) { if (SetBits (nTaxa) == ERROR) return (ERROR); } /* get size of unsigned */ unsignedIntsize = sizeof(unsigned); /* get down pass sequence */ i = 0; GetUserDownPass (sumTreeRoot, sumTreeDP, &i); /* pass down tree, assigning bits for each node */ for (i=0; iindex)*nUnsignedNeeded + j] = 0; if (p->left == NULL && p->right == NULL && p->anc != NULL) { offSet = 0; while ((p->index+1) > unsignedIntsize*8*(offSet+1)) offSet++; x = 1 << (p->index - offSet*unsignedIntsize*8); treeBits[(p->index)*nUnsignedNeeded + offSet] = x; } else if (p->left != NULL && p->right == NULL && p->anc == NULL && isSumTreeRooted == NO) { offSet = 0; while ((p->index+1) > unsignedIntsize*8*(offSet+1)) offSet++; x = 1; x <<= (p->index - offSet*unsignedIntsize*8); treeBits[(p->index)*nUnsignedNeeded + offSet] = x; } else if (p->left != NULL && p->right != NULL && p->anc != NULL) { for (j=0; jindex)*nUnsignedNeeded + j] = treeBits[(p->left->index)*nUnsignedNeeded + j] | treeBits[(p->right->index)*nUnsignedNeeded + j]; } } /*printf ("%4d: ", p->index); for (j=nUnsignedNeeded-1; j>=0; j--) DisplayBits (treeBits[(p->index)*nUnsignedNeeded + j]); printf ("\n");*/ } /* flip bits, if necessary (first taxon is always 0) */ if (isSumTreeRooted == NO) { for (i=0; iindex); for (j=nUnsignedNeeded-1; j>=0; j--) DisplayBits (treeBits[(p->index)*nUnsignedNeeded + j]); printf ("\n");*/ } } } } /* tally each partition and branch length */ if (numSumTrees > burnIn) { numSumTreesSampled++; for (i=0; ianc != NULL) { if (p->anc->anc != NULL) { foundPartition = NO; for (k=0; kindex)*nUnsignedNeeded + j] != treePartsFound[k*nUnsignedNeeded + j]) { allSame = NO; break; } } if (allSame == YES) { numFoundOfThisPart[k]++; foundPartition = YES; aOldBrlens[k] = aBrlens[k]; aBrlens[k] = aOldBrlens[k] + (p->length - aOldBrlens[k]) / (numFoundOfThisPart[k]); sBrlens[k] = sBrlens[k] + (p->length - aOldBrlens[k]) * (p->length - aOldBrlens[k]); break; } } if (foundPartition == NO) { if (numTreePartsFound+1 == numPartsAllocated) { if (RealloateTreeParts () == ERROR) { printf (" ERROR: Could not increase size of partitions.\n"); return (ERROR); } } for (j=0; jindex)*nUnsignedNeeded + j]; numFoundOfThisPart[numTreePartsFound] = 1; aBrlens[numTreePartsFound] = p->length; sBrlens[numTreePartsFound] = 0.0; numTreePartsFound++; } } } else { if (isSumTreeRooted == NO) { foundPartition = NO; for (k=0; kindex)*nUnsignedNeeded + j] != treePartsFound[k*nUnsignedNeeded + j]) { allSame = NO; break; } } if (allSame == YES) { numFoundOfThisPart[k]++; foundPartition = YES; aOldBrlens[k] = aBrlens[k]; aBrlens[k] = aOldBrlens[k] + (p->left->length - aOldBrlens[k]) / (numFoundOfThisPart[k]); sBrlens[k] = sBrlens[k] + (p->left->length - aOldBrlens[k]) * (p->left->length - aOldBrlens[k]); break; } } if (foundPartition == NO) { if (numTreePartsFound+1 == numPartsAllocated) { if (RealloateTreeParts () == ERROR) { printf (" ERROR: Could not increase size of partitions.\n"); return (ERROR); } } for (j=0; jindex)*nUnsignedNeeded + j]; numFoundOfThisPart[numTreePartsFound] = 1; aBrlens[numTreePartsFound] = p->left->length; sBrlens[numTreePartsFound] = 0.0; numTreePartsFound++; } } } } } # if 0 printf (" Tree %d:\n", numSumTrees); if (ShowTree (sumTreeRoot, isSumTreeRooted) == ERROR) { return (ERROR); } # endif # if 0 ShowNodes (sumTreeRoot, 2, isSumTreeRooted); # endif break; case CMD_SUM_TREE8: /* */ sscanf (token, "%lf", &tempD); if (pPtr->marked == NO) pPtr->length = tempD; else { if (pPtr->left != NULL) pPtr->left->length = tempD; else { printf (" ERROR: Cannot assign branch length to left node\n"); printf ("token = %s\n", token); return (ERROR); } } foundBrlens = YES; break; case CMD_INIT_TREE_BLK: isTranslated = NO; numSumTrees = 0; numSumTreesSampled = 0; break; case CMD_SUM_TREE_END: printf (" Read %d trees from file, %d of which were sampled\n", numSumTrees, numSumTreesSampled); if (numSumTreesSampled == 0) { printf ("\n ERROR: No trees read in.\n\n", conFile); free (sumTreeDP); free (treeBits); free (treePartsFound); free (numFoundOfThisPart); free (aBrlens); free (aOldBrlens); free (sBrlens); free (taxonMask); ErrorExit (); return (ERROR); } unsignedIntsize = sizeof(unsigned); /* sort partitions */ SortParts (numFoundOfThisPart, numTreePartsFound); /* print partitions */ strcpy (conFile, sumTreeFile); strcat (conFile, ".parts"); if ((fp = fopen (conFile, "w")) == NULL) { printf ("\n ERROR: Could not open file \"%s\"\n\n", conFile); free (sumTreeDP); free (treeBits); free (treePartsFound); free (numFoundOfThisPart); free (aBrlens); free (aOldBrlens); free (sBrlens); free (taxonMask); ErrorExit (); return (ERROR); } printf ("\n"); printf (" General explanation:\n\n"); printf (" A taxon bipartition is specified by removing a\n"); printf (" branch, dividing the species into those to the\n"); printf (" left and those to the right of the branch. Here,\n"); printf (" taxa to one side of the removed branch are denoted\n"); printf (" \".\" and those to the other side are denoted \"*\".\n"); printf (" The output includes the bipartition number (sorted\n"); printf (" from highest to lowest probability), bipartition\n"); printf (" (e.g., ...**..), number of times the bipartition\n"); printf (" was observed, the posterior probability of the \n"); printf (" bipartition, and, if branch lengths were recorded \n"); printf (" on the trees in the file, the average and variance\n"); printf (" of the branch lengths. Each \".\" or \"*\" in the\n"); printf (" bipartition represents a taxon that is to the left \n"); printf (" or right of the removed branch. A list of the taxa \n"); printf (" in the bipartitions is listed before the list of \n"); printf (" bipartitions. \n\n"); printf (" List of taxa in bipartitions:\n\n"); for (k=0; k= unsignedIntsize*8*(offSet+1)) { offSet++; x = 1; } } if (foundBrlens == YES) { printf (" -- %5d %1.6lf %1.8lf (%1.8lf)\n", numFoundOfThisPart[k], (double)numFoundOfThisPart[k]/numSumTreesSampled, aBrlens[k], (double)sBrlens[k]/(numSumTreesSampled-1)); fprintf (fp, "\t%d\t%1.6lf\t%1.8lf\t%1.8lf\n", numFoundOfThisPart[k], (double)numFoundOfThisPart[k]/numSumTreesSampled, aBrlens[k], (double)sBrlens[k]/(numSumTreesSampled-1)); } else { printf (" -- %5d %1.6lf\n", numFoundOfThisPart[k], (double)numFoundOfThisPart[k]/numSumTreesSampled); fprintf (fp, "\t%d\t%1.6lf\n", numFoundOfThisPart[k], (double)numFoundOfThisPart[k]/numSumTreesSampled); } } # if !defined (CONTREE) fclose (fp); # endif # if defined (CONTREE) /* allocate space for consensus tree nodes */ conNodes = (PolyNode *)malloc((size_t) (2 * nTaxa * sizeof(PolyNode))); if (!conNodes) { printf (" ERROR: Could not allocate conNodes\n"); free (sumTreeDP); free (treeBits); free (treePartsFound); free (numFoundOfThisPart); free (aBrlens); free (aOldBrlens); free (sBrlens); free (taxonMask); fclose(fp); ErrorExit (); return (ERROR); } /* transform taxonLabels */ for (i=0; ix counts number of subtended terminals */ /* make sure conRoot->left is in outgroup */ conRoot = &conNodes[nTaxa]; conRoot->anc = conRoot->sib = NULL; conRoot->x = nTaxa; j = FirstTaxonInPartition(outgroupPartition, nUnsignedNeeded); 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 more frequent to less frequent */ /* for quit test to work when a 50% maj rule tree is requested */ nextConNode = nTaxa + 1; if (isSumTreeRooted) targetNode = 2*nTaxa - 2; else targetNode = 2*nTaxa - 3; for (i=0; i targetNode) break; freq = (double)numFoundOfThisPart[i]/ (double)numSumTreesSampled; if (freq < 0.50 && consensusType == HALF_COMPAT) break; // get partition partition = &treePartsFound[i*nUnsignedNeeded]; // flip bits if necessary if(!isSumTreeRooted) { if (!IsPartNested(partition, ingroupPartition, nUnsignedNeeded) && !IsPartCompatible(partition, ingroupPartition, nUnsignedNeeded)) FlipBits(partition, nUnsignedNeeded, lastMask); } // count bits in this partition for (j=nBits=0; j 1) // this is an informative partition { // find anc of partition j = FirstTaxonInPartition (partition, nUnsignedNeeded); 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, nUnsignedNeeded)) isCompat = NO; for (q=cp->left; q!=NULL; q=q->sib) { if (q->x > 1 && !IsPartCompatible(q->partition, partition, nUnsignedNeeded)) isCompat = NO; if (isCompat == NO) break; } if (isCompat == NO) continue; // set new node q = &conNodes[nextConNode++]; q->length = aBrlens[i]; 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, nUnsignedNeeded)) || (r->x == 1 && (partition[r->index / (unsignedIntsize * 8)] & (1 << (r->index % (unsignedIntsize * 8)))) != 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, nUnsignedNeeded); q = &conNodes[j]; q->length = aBrlens[i]; q->support = freq * 100; } } /* draw tree to stdout and fp */ fprintf(fp,"\nClade credibility values:\n\n"); printf("\nClade credibility values:\n\n"); ShowConTree (stdout, nextConNode, conRoot, 80); ShowConTree (fp, nextConNode, conRoot, 80); if (foundBrlens == YES) { fprintf(fp,"Phylogram:\n\n"); printf("Phylogram:\n\n"); ShowConPhylogram (stdout, nextConNode, conRoot, 80); ShowConPhylogram (fp, nextConNode, conRoot, 80); } fclose (fp); /* reset taxonLabels */ for (i=0; i= 2*nTaxa) { printf (" ERROR: Too many compatible partitions found.\n"); free (constraintTree); free (taxaNotInConstraint); free (taxaInConstraint); free (compatParts); free (tempParts); free (numOnes); free (partLengths); free (sumTreeDP); free (treeBits); free (treePartsFound); free (numFoundOfThisPart); free (aBrlens); free (aOldBrlens); free (sBrlens); free (taxonMask); ErrorExit (); return (ERROR); } for (j=0; j= unsignedIntsize*8*(offSet+1)) { offSet++; x = 1; } } numOnes[i] = k; } /* sort partitions from largest to smallest */ SortCompatParts (numOnes, compatParts, partLengths, nCompatParts); # if 0 for (k=0; k= unsignedIntsize*8*(offSet+1)) { offSet++; x = 1; } } printf (" -- %5d %lf\n", numOnes[k], partLengths[k]); } # endif /* set up initial star tree */ conRoot = &constraintTree[2*nTaxa-1]; cq = NULL; for (i=0; ilength = 0.0; cp->anc = conRoot; if (cq != NULL) cq->rightSib = cp; else cp->anc->left = cp; cq = cp; } nextConNode = nTaxa; /* build up constraint tree with compatible partitions */ for (i=0; i= unsignedIntsize*8*(offSet+1)) { offSet++; x = 1; } } } # if 0 printf ("in: "); for (j=0; jindex); printf ("\n"); # endif /* get pointer to ancestor */ cp = taxaInConstraint[0]->anc; isCompat = YES; for (j=0; janc != cp) { isCompat = NO; break; } } if (isCompat == NO) { printf (" ERROR: Strange incompatibility in consensus tree found.\n"); free (constraintTree); free (taxaNotInConstraint); free (taxaInConstraint); free (compatParts); free (tempParts); free (numOnes); free (partLengths); free (sumTreeDP); free (treeBits); free (treePartsFound); free (numFoundOfThisPart); free (aBrlens); free (aOldBrlens); free (sBrlens); free (taxonMask); ErrorExit (); return (ERROR); } //printf ("cp = %d (", cp->index); /* get old left-most node */ cl = cp->left; /* get old right-most node */ for (cq = cp->left; cq != NULL; cq = cq->rightSib) { if (cq != NULL) if (cq->rightSib == NULL) cr = cq; //printf ("%d ", cq->index); } //printf (")\n"); //printf ("cr = %d\n", cr->index); /* make a list of nodes not in constraint */ nTaxaNotInConstraint = 0; for (cq = cp->left; cq != NULL; cq = cq->rightSib) { tempInt = NO; for (j=0; jindex); printf ("\n"); # endif /* set up constraint */ if (numOnes[i] > 1) { /* get new node */ if (nextConNode + 1 < 2*nTaxa - 1) { cn = &constraintTree[nextConNode++]; cn->length = partLengths[i]; } else { printf ("new = PROBLEM!\n"); } cp->left = cn; cn->anc = cp; if (nTaxaNotInConstraint == 0) cn->rightSib = NULL; else { cn->rightSib = taxaNotInConstraint[0]; for (j=0; jrightSib = taxaNotInConstraint[j+1]; else cq->rightSib = NULL; } } cn->left = taxaInConstraint[0]; for (j=0; janc = cn; if (j+1 < nTaxaInConstraint) cq->rightSib = taxaInConstraint[j+1]; else cq->rightSib = NULL; } } else { offSet = 0; x = 1; for (j=0; jlength = partLengths[i]; break; } x <<= 1; if (j+1 >= unsignedIntsize*8*(offSet+1)) { offSet++; x = 1; } } } } /*for (i=0; i<2*nTaxa; i++) { cp = &constraintTree[i]; printf ("%4d (%d %d %d)\n", Dex3(cp), Dex3(cp->left), Dex3(cp->rightSib), Dex3(cp->anc)); } printf ("conRoot = %d\n", conRoot->index);*/ #endif /* write constraint tree */ strcpy (conFile, sumTreeFile); strcat (conFile, ".con"); if ((fp = fopen (conFile, "w")) == NULL) { printf ("\n ERROR: Could not open file \"%s\"\n\n", conFile); # if !defined (CONTREE) free (constraintTree); free (taxaNotInConstraint); free (taxaInConstraint); free (compatParts); free (tempParts); free (numOnes); free (partLengths); # else free (conNodes); free (outgroupPartition); # endif free (sumTreeDP); free (treeBits); free (treePartsFound); free (numFoundOfThisPart); free (aBrlens); free (aOldBrlens); free (sBrlens); free (taxonMask); ErrorExit (); return (ERROR); } fprintf (fp, "#NEXUS\n\n"); fprintf (fp, "begin trees;\n"); if (isTranslated == YES) { fprintf (fp, "translate\n"); for (j=0; j maxX) maxX = paramInfo[i * nColumns + whichIsX]; if (paramInfo[i * nColumns + whichIsY] > maxY) maxY = paramInfo[i * nColumns + whichIsY]; } //printf ("%lf %lf %lf %lf\n", minX, maxX, minY, maxY); for (i=0; i= SCREENWIDTH) k = SCREENWIDTH - 1; meanY[k] += paramInfo[i * nColumns + whichIsY]; numY[k]++; } printf ("\n\n +"); for (i=0; i=0; j--) { printf (" |"); 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) printf ("*"); else printf (" "); } else { printf (" "); } } printf ("|\n"); } printf (" +"); for (i=0; i= paramBurn) tempInfo[k++] = paramInfo[i * nColumns + j]; } /* sort the contents of the vector */ SortVector (tempInfo, k); /*for (i=0; i= HEADER_SIZE) { printf (" ERROR: The header is too large\n"); free (headerNames); ErrorExit (); } headerNames[headerSize++] = token[i]; i++; } nColumns++; break; case CMD_SUM_PARM4: if (headerSize >= HEADER_SIZE) { printf (" ERROR: The header is too large\n"); free (headerNames); ErrorExit (); } headerNames[headerSize++] = '|'; break; case CMD_SUM_PARM5: if (headerSize+2 >= HEADER_SIZE) { printf (" ERROR: The header is too large\n"); free (headerNames); ErrorExit (); } headerNames[headerSize++] = '|'; headerNames[headerSize++] = '\0'; //printf ("header = %s\n", headerNames); /* allocate memory for data */ paramInfo = (double *)malloc((size_t) (100 * nColumns * sizeof(double))); if (!paramInfo) { printf (" ERROR: Could not allocate paramInfo\n"); free (headerNames); ErrorExit (); return (ERROR); } curNumRows = 100; //printf ("number of columns = %d\n", nColumns); whichColumn = whichRow = 0; break; case CMD_SUM_PARM6: /* data */ tempNumColumns = 0; break; case CMD_SUM_PARM7: sscanf (token, "%lf", &tempD); if (whichRow == curNumRows) { /* reallocate paraminfo */ curNumRows += 100; paramInfo = (double *)realloc(paramInfo, (size_t) (curNumRows * nColumns * sizeof(double))); if (!paramInfo) { printf (" ERROR: Could not allocate paramInfo\n"); free (headerNames); free (paramInfo); ErrorExit (); return (ERROR); } } if (whichColumn == nColumns) { printf (" ERROR: Too many columns\n"); free (headerNames); free (paramInfo); ErrorExit (); return (ERROR); } paramInfo[whichRow * nColumns + whichColumn] = tempD; whichColumn++; break; case CMD_SUM_PARM8: break; case CMD_SUM_PARM9: whichColumn = 0; whichRow++; break; case CMD_SUM_PARM10: sscanf (token, "%lf", &tempD); if (whichRow == curNumRows) { /* reallocate paraminfo */ curNumRows += 100; paramInfo = (double *)realloc(paramInfo, (size_t) (curNumRows * nColumns * sizeof(double))); if (!paramInfo) { printf (" ERROR: Could not allocate paramInfo\n"); free (headerNames); free (paramInfo); ErrorExit (); return (ERROR); } } if (whichColumn == nColumns) { printf (" ERROR: Too many columns\n"); free (headerNames); free (paramInfo); ErrorExit (); return (ERROR); } paramInfo[whichRow * nColumns + whichColumn] = tempD * (-1.0); whichColumn++; break; case CMD_AUTOCLOSE1: autoclose = YES; printf (" Setting autoclose=yes\n"); break; case CMD_AUTOCLOSE2: autoclose = NO; printf (" Setting autoclose=no\n"); break; case CMD_NPERB: sscanf (token, "%d", &tempInt); nPerturbations = tempInt; break; case CMD_MCMCP: printf (" Setting chain parameters\n"); break; case CMD_SET_AAMODEL1: aaModelType = POISSON; printf (" Setting aamodel=poisson\n"); break; case CMD_SET_AAMODEL2: aaModelType = EQUALIN; printf (" Setting aamodel=equalin\n"); break; case CMD_SET_AAMODEL3: aaModelType = JONES; printf (" Setting aamodel=jones\n"); break; case CMD_SET_AAMODEL4: aaModelType = DAYHOFF; printf (" Setting aamodel=dayhoff\n"); break; case CMD_SET_AAMODEL5: aaModelType = GTR; printf (" Setting aamodel=gtr\n"); break; case CMD_SET_ANCFILE1: saveAncFiles = YES; printf (" Setting ancfile=yes\n"); break; case CMD_SET_ANCFILE2: saveAncFiles = NO; printf (" Setting ancfile=no\n"); break; case CMD_SET_CODON1: enforceCodonModel = YES; printf (" Setting enforcecodon=yes\n"); break; case CMD_SET_CODON2: enforceCodonModel = NO; printf (" Setting enforcecodon=no\n"); break; case CMD_SET_OMEGAVAR1: strcpy (codonModelType, token); printf (" Setting omegavar=equal\n"); break; case CMD_SET_OMEGAVAR2: strcpy (codonModelType, token); printf (" Setting omegavar=ny98\n"); break; case CMD_SET_OMEGAVAR3: strcpy (codonModelType, token); printf (" Setting omegavar=ac1ny98\n"); break; case CMD_SET_OMEGAVAR4: strcpy (codonModelType, token); printf (" Setting omegavar=ac2ny98\n"); break; case CMD_SET_OMEGAVAR5: strcpy (codonModelType, token); printf (" Setting omegavar=covny98\n"); break; case CMD_SET_OMEGA: if (token[0] == 'e') { strcpy(omegaRatioModel, token); omegaRatio = 1.0; } else { strcpy (omegaRatioModel, "user defined"); sscanf (token, "%lf", &tempD); omegaRatio = tempD; } printf (" Setting omega = %lf\n", omegaRatio); break; case CMD_OMEGAPR1: strcpy (omegaprModel, "fixed"); printf (" Fixing omega\n"); break; case CMD_OMEGAPR2: strcpy (omegaprModel, "uniform"); printf (" Uniform prior on omega\n"); break; case CMD_OMEGAPR3: sscanf (token, "%lf", &tempD); if (tempD < 0.0) { printf (" ERROR: Lower bound on uniform must be positive\n"); return (ERROR); } omegaprUni[0] = tempD; break; case CMD_OMEGAPR4: sscanf (token, "%lf", &tempD); if (tempD <= omegaprUni[0]) { printf (" ERROR: Upper bound on uniform must be greater than lower bound\n"); return (ERROR); } omegaprUni[1] = tempD; break; case CMD_OMEGAPR5: strcpy (omegaprModel, "exponential"); printf (" Exponential prior on omega\n"); break; case CMD_OMEGAPR6: sscanf (token, "%lf", &tempD); if (tempD <= 0.0) { printf (" ERROR: Parameter of exponential must be greater than 0.\n"); return (ERROR); } omegaprExp = tempD; break; case CMD_SET_CODE: if (token[0] == 'u') aaCode = UNIV; else if (token[0] == 'v') aaCode = VERTMT; else if (token[0] == 'm' && token[1] == 'y') aaCode = MYCO; else if (token[0] == 'y') aaCode = YEAST; else if (token[0] == 'c') aaCode = CILIATES; else aaCode = METMT; printf (" Setting code to %s\n", token); break; case CMD_INFERSEL_YES: printf (" Setting inferpossel=yes\n"); inferPoseSel = YES; break; case CMD_INFERSEL_NO: printf (" Setting inferpossel=no\n"); inferPoseSel = NO; break; case CMD_INIT_LAG1: strcpy (lagModel, "user defined"); sscanf (token, "%d", &tempInt); lagNumber = tempInt; printf (" Setting user-defined lag of %d\n", lagNumber); break; case CMD_INIT_LAG2: strcpy (lagModel, token); printf (" Estimating lag\n"); break; case CMD_SET_LAG_PR1: strcpy (lagprModel, "fixed"); printf (" Fixing lag\n"); strcpy (lagprModel, "fixed"); break; case CMD_SET_LAG_PR2: sscanf (token, "%lf", &tempD); lagprUni[0] = tempD; break; case CMD_SET_LAG_PR3: sscanf (token, "%lf", &tempD); if (tempD <= lagprUni[0]) { printf (" ERROR: Upper bound on uniform must be greater than lower bound\n"); return (ERROR); } lagprUni[1] = tempD; printf (" Uniform (%1.3lf, %1.3lf) prior on lag\n", lagprUni[0], lagprUni[1]); strcpy (lagprModel, "uniform"); break; case CMD_SET_LAG_PR4: sscanf (token, "%lf", &tempD); if (tempD <= 0.0) { printf (" ERROR: Exponential parameter must be greater than 0\n"); return (ERROR); } lagprExp = tempD; printf (" Exponential (%1.3lf) prior on lag\n", lagprExp); strcpy (lagprModel, "exponential"); break; case CMD_SET_INFERRATES1: printf (" Setting inferrates=yes\n"); inferRates = YES; break; case CMD_SET_INFERRATES2: printf (" Setting inferrates=no\n"); inferRates = NO; break; case CMD_SET_COVARION1: printf (" Setting covarion=yes\n"); useCovarion = YES; break; case CMD_SET_COVARION2: printf (" Setting covarion=no\n"); useCovarion = NO; break; case CMD_SET_PROPS1: printf (" Setting %s\n", token); if (!strcmp("qmatwin", token)) propSet = 0; else if (!strcmp("pialpha", token)) propSet = 1; else if (!strcmp("ratewin", token)) propSet = 2; else if (!strcmp("tuning", token)) propSet = 3; else if (!strcmp("clocktune", token)) propSet = 4; else if (!strcmp("timetune", token)) propSet = 5; else if (!strcmp("mtune", token)) propSet = 6; else if (!strcmp("extendp", token)) propSet = 7; else if (!strcmp("wormtune", token)) propSet = 8; else if (!strcmp("bdwin", token)) propSet = 9; else if (!strcmp("omegawin", token)) propSet = 10; else if (!strcmp("rhowin", token)) propSet = 11; else if (!strcmp("omegaalpha", token)) propSet = 12; else if (!strcmp("invslide", token)) propSet = 13; else if (!strcmp("switchwin", token)) propSet = 14; else if (!strcmp("alphawin", token)) propSet = 15; else if (!strcmp("reconlim", token)) propSet = 16; else if (!strcmp("tbrextendp", token)) propSet = 17; else if (!strcmp("tbrtune", token)) propSet = 18; else if (!strcmp("bltune", token)) propSet = 19; break; case CMD_SET_PROPS2: sscanf (token, "%lf", &tempD); if (propSet == 0) qMatWinProp = tempD; else if (propSet == 1) piAlphaProp = tempD; else if (propSet == 2) rateWinProp = tempD; else if (propSet == 3) tuningProp = tempD; else if (propSet == 4) clockTuneProp = tempD; else if (propSet == 5) timeTuneProp = tempD; else if (propSet == 6) mTuneProp = tempD; else if (propSet == 7) extendPProp = tempD; else if (propSet == 8) wormTuneProp = tempD; else if (propSet == 9) bdWinProp = tempD; else if (propSet == 10) omegaWinProp = tempD; else if (propSet == 11) rhoWinProp = tempD; else if (propSet == 12) omegaAlphaProp = tempD; else if (propSet == 13) invSlideProp = tempD; else if (propSet == 14) switchWinProp = tempD; else if (propSet == 15) alphaWinProp = tempD; else if (propSet == 16) reconLimProp = tempD; else if (propSet == 17) tbrExtendPProp = tempD; else if (propSet == 18) tbrTuneProp = tempD; else if (propSet == 19) blTuneProp = tempD; break; case CMD_SET_COVTYPE: strcpy (covarionModelType, token); printf (" Setting covtype=%s\n", covarionModelType); break; case CMD_SET_CODING: strcpy (codingType, token); printf (" Setting coding=%s\n", codingType); break; case CMD_SET_STATE_PR1: strcpy (stateFreqPrModel, "fixed"); printf (" State frequencies fixed\n"); break; case CMD_SET_STATE_PR2: sscanf (token, "%lf", &tempD); statefreqprDir = tempD; strcpy (stateFreqPrModel, "dirichlet"); printf (" State frequencies follow a dirichlet(%lf)\n", tempD); break; case CMD_SET_MOVE_PROPS: printf (" Setting move rates\n\n"); for (;;) { printf (" 1: %1.2lf = relative probability of changes to Q matrix\n", moveRates[0]); printf (" 2: %1.2lf = relative probability of changes to base frequencies\n", moveRates[1]); printf (" 3: %1.2lf = relative probability of changes to gamma shape\n", moveRates[2]); printf (" 4: %1.2lf = relative probability of changes to site rates\n", moveRates[3]); printf (" 5: %1.2lf = relative probability of changes to tree (stochastic NNI)\n", moveRates[4]); printf (" 6: %1.2lf = relative probability of changes to clock tree\n", moveRates[5]); printf (" 7: %1.2lf = relative probability of changes to time constrained tree\n", moveRates[6]); printf (" 8: %1.2lf = relative probability of changes to tree height\n", moveRates[7]); printf (" 9: %1.2lf = relative probability of changes to tree (one worm)\n", moveRates[8]); printf (" 10: %1.2lf = relative probability of changes to node time\n", moveRates[9]); printf (" 11: %1.2lf = relative probability of changes to bd parameters\n", moveRates[10]); printf (" 12: %1.2lf = relative probability of changes to omega\n", moveRates[11]); printf (" 13: %1.2lf = relative probability of changes to omega probabilities\n", moveRates[12]); printf (" 14: %1.2lf = relative probability of changes to auto correlation\n", moveRates[13]); printf (" 15: %1.2lf = relative probability of changes to lag\n", moveRates[14]); printf (" 16: %1.2lf = relative probability of changes to tree (stochastic SPR)\n", moveRates[15]); printf (" 17: %1.2lf = relative probability of changes to proportion invariant\n", moveRates[16]); printf (" 18: %1.2lf = relative probability of changes to covarion parameter\n", moveRates[17]); printf (" 19: %1.2lf = relative probability of hanges to tree (TBR1)\n", moveRates[18]); printf (" 20: %1.2lf = relative probability of changes to GC swith parameters\n", moveRates[19]); printf (" 21: %1.2lf = relative probability of changes to tree (TBR2)\n", moveRates[20]); printf (" 22: %1.2lf = relative probability of changes to a branch length\n", moveRates[21]); printf (" 23: %1.2lf = relative probability of changes to tree (eraser)\n", moveRates[22]); printf (" 24: %1.2lf = relative probability of changes to tree (biased SPR)\n", moveRates[23]); printf (" 25: %1.2lf = relative probability of changes to tree (biased TBR)\n", moveRates[24]); printf ("\n"); printf (" Which move do you wish to change (0 to exit): "); tempInt = (int)EnterNum (); if (tempInt == 0) break; else if (tempInt < 0 || tempInt > 25) { printf (" You must enter a number between 0 and 25\n\n"); } else { printf (" Enter new relative proposal probability: "); tempD = EnterNum (); moveRates[tempInt - 1] = tempD; } } break; case CMD_SET_MCMC_BURN: sscanf (token, "%d", &tempInt); mcmcBurn = tempInt; break; case CMD_SET_CYCLEFREQ: sscanf (token, "%d", &tempInt); cycleFreq = tempInt; break; case CMD_SET_CALCPSEUDOBF: if (!strcmp("yes", token)) calcPseudoBf = YES; else calcPseudoBf = NO; break; case CMD_SET_PARS_CRIT: printf (" Setting parsmodel=%s\n", token); if (!strcmp("yes", token)) useParsCriterion = YES; else useParsCriterion = NO; break; case CMD_SET_TYPESET1: //printf ("typset1 (%s)\n", token); if (isMatrixAllocated == NO) { printf (" ERROR: Data not read in\n"); return (ERROR); } break; case CMD_SET_TYPESET2: //printf ("typset2 (%s)\n", token); i = 0; while (token[i] != '\0') { typeDefLabels[nTypeDefs][i] = token[i]; i++; } typeDefLabels[nTypeDefs][i] = '|'; charSetFrom = -1; charSetTo = -1; charSetWait = NO; charSetSlashWait = NO; charSetPotWait = YES; charSetFreq = 1; printf (" Assigning typeset set \""); i = 0; while (typeDefLabels[nTypeDefs][i] != '|') { printf ("%c", typeDefLabels[nTypeDefs][i]); i++; } printf ("\"\n"); if (nTypeDefs > 0) { typeDef = (int *)realloc(typeDef, (size_t) ((nTypeDefs+1) * nChar * sizeof(int))); if (!typeDef) { printf (" ERROR: Could not reallocate typeDef\n"); return (ERROR); } for (i=nTypeDefs*nChar; i<(nTypeDefs+1)*nChar; i++) typeDef[i] = UNORDERED; } break; case CMD_SET_TYPESET3: //printf ("typset3 (%s)\n", token); printf (" ERROR: You cannot use typeset names that are digits\n"); return (ERROR); break; case CMD_SET_TYPESET4: //printf ("typset4 (%s)\n", token); if (!strcmp("ord", token)) charOrdering = ORDERED; else charOrdering = UNORDERED; break; case CMD_SET_TYPESET5: //printf ("typset5 (%s)\n", token); break; case CMD_SET_TYPESET6: //printf ("typset6 (%s)\n", token); if (charSetWait == NO && charSetSlashWait == NO && charSetPotWait == NO) { j = charSetFreq; for (i=charSetFrom-1; i<=charSetTo-1; i++) { if (j % charSetFreq == 0) typeDef[nTypeDefs*nChar + i] = charOrdering; j++; } charSetFrom = -1; charSetTo = -1; charSetWait = NO; charSetSlashWait = NO; charSetFreq = 1; charSetPotWait = YES; } sscanf (token, "%d", &tempInt); if (tempInt <= 0 || tempInt > nChar) { printf (" ERROR: Number (%d) out of range\n", tempInt); return (ERROR); } if (charSetSlashWait == YES && charSetWait == YES) { printf (" ERROR: Problem assigning character set\n"); return (ERROR); } if (charSetSlashWait == YES) { charSetFreq = tempInt; charSetSlashWait = NO; charSetPotWait = NO; } else if (charSetWait == YES) { charSetTo = tempInt; charSetWait = NO; charSetPotWait = NO; } else { if (charSetWait == NO) { charSetFrom = tempInt; charSetTo = tempInt; charSetPotWait = NO; } else { charSetTo = tempInt; charSetPotWait = NO; } } break; case CMD_SET_TYPESET7: //printf ("typset7 (%s)\n", token); charSetWait = YES; charSetPotWait = YES; break; case CMD_SET_TYPESET8: //printf ("typset8 (%s)\n", token); charSetSlashWait = YES; charSetPotWait = YES; break; case CMD_SET_TYPESET9: //printf ("typset9 (%s)\n", token); if (charSetWait == NO && charSetSlashWait == NO && charSetPotWait == NO) { j = charSetFreq; for (i=charSetFrom-1; i<=charSetTo-1; i++) { if (j % charSetFreq == 0) typeDef[nTypeDefs*nChar + i] = charOrdering; j++; } charSetFrom = -1; charSetTo = -1; charSetWait = NO; charSetSlashWait = NO; charSetFreq = 1; charSetPotWait = YES; } /*for (i=0; i<(nTypeDefs+1)*nChar; i++) { printf ("%d", typeDef[i]); } printf ("\n");*/ nTypeDefs++; break; case CMD_SET_TYPESET10: //printf ("typset10 (%s)\n", token); foundSet = NO; for (i=0; i 0) { nStatesDef = (int *)realloc(nStatesDef, (size_t) ((nNStatesDefs+1) * nChar * sizeof(int))); if (!nStatesDef) { printf (" ERROR: Could not reallocate nStatesDef\n"); return (ERROR); } for (i=nNStatesDefs*nChar; i<(nNStatesDefs+1)*nChar; i++) nStatesDef[i] = nCharStates; } break; case CMD_SET_STATESSET3: //printf ("typset3 (%s)\n", token); printf (" ERROR: You cannot use nstatesset names that are digits\n"); return (ERROR); break; case CMD_SET_STATESSET4: //printf ("typset4 (%s)\n", token); sscanf (token, "%d", &charOrdering); break; case CMD_SET_STATESSET5: //printf ("typset5 (%s)\n", token); break; case CMD_SET_STATESSET6: //printf ("typset6 (%s)\n", token); if (charSetWait == NO && charSetSlashWait == NO && charSetPotWait == NO) { j = charSetFreq; for (i=charSetFrom-1; i<=charSetTo-1; i++) { if (j % charSetFreq == 0) nStatesDef[nNStatesDefs*nChar + i] = charOrdering; j++; } charSetFrom = -1; charSetTo = -1; charSetWait = NO; charSetSlashWait = NO; charSetFreq = 1; charSetPotWait = YES; } sscanf (token, "%d", &tempInt); if (tempInt <= 0 || tempInt > nChar) { printf (" ERROR: Number (%d) out of range\n", tempInt); return (ERROR); } if (charSetSlashWait == YES && charSetWait == YES) { printf (" ERROR: Problem assigning character set\n"); return (ERROR); } if (charSetSlashWait == YES) { charSetFreq = tempInt; charSetSlashWait = NO; charSetPotWait = NO; } else if (charSetWait == YES) { charSetTo = tempInt; charSetWait = NO; charSetPotWait = NO; } else { if (charSetWait == NO) { charSetFrom = tempInt; charSetTo = tempInt; charSetPotWait = NO; } else { charSetTo = tempInt; charSetPotWait = NO; } } break; case CMD_SET_STATESSET7: //printf ("typset7 (%s)\n", token); charSetWait = YES; charSetPotWait = YES; break; case CMD_SET_STATESSET8: //printf ("typset8 (%s)\n", token); charSetSlashWait = YES; charSetPotWait = YES; break; case CMD_SET_STATESSET9: //printf ("typset9 (%s)\n", token); if (charSetWait == NO && charSetSlashWait == NO && charSetPotWait == NO) { j = charSetFreq; for (i=charSetFrom-1; i<=charSetTo-1; i++) { if (j % charSetFreq == 0) nStatesDef[nNStatesDefs*nChar + i] = charOrdering; j++; } charSetFrom = -1; charSetTo = -1; charSetWait = NO; charSetSlashWait = NO; charSetFreq = 1; charSetPotWait = YES; } /*for (i=0; i<(nNStatesDefs+1)*nChar; i++) { printf ("%d", nStatesDef[i]); } printf ("\n");*/ nNStatesDefs++; break; case CMD_SET_STATESSET10: //printf ("typset10 (%s)\n", token); foundSet = NO; for (i=0; ileft, i, isThisTreeRooted); FinishTree (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)++; } } } #if defined (CONTREE) int FirstTaxonInPartition (unsigned int *partition, int length) { int i, j, nBits, taxon; unsigned int x; nBits = sizeof(unsigned int) * 8; for (i=taxon=0; ileft != NULL) { for (q=p->left; q!=NULL; q=q->sib) GetConDownPass(downPass, q, i); } downPass[(*i)++] = p; } #endif void GetToken (int *tokenType) { int allNumbers; register char *temp; (*tokenType) = 0; temp = token; while (IsWhite(*tokenP) == 1 || IsWhite(*tokenP) == 2) { if (IsWhite(*tokenP) == 2) *tokenType = RETURNSYMBOL; ++tokenP; } if (IsIn(*tokenP,"=")) { *temp++ = *tokenP++; *tokenType = EQUALSIGN; } else if (IsIn(*tokenP,";")) { *temp++ = *tokenP++; *tokenType = SEMICOLON; } else if (IsIn(*tokenP,":")) { *temp++ = *tokenP++; *tokenType = COLON; } else if (IsIn(*tokenP,",")) { *temp++ = *tokenP++; *tokenType = COMMA; } else if (IsIn(*tokenP,"#")) { *temp++ = *tokenP++; *tokenType = POUNDSIGN; } else if (IsIn(*tokenP,"(")) { *temp++ = *tokenP++; *tokenType = LEFTPAR; } else if (IsIn(*tokenP,")")) { *temp++ = *tokenP++; *tokenType = RIGHTPAR; } else if (IsIn(*tokenP,"[")) { *temp++ = *tokenP++; *tokenType = LEFTCOMMENT; } else if (IsIn(*tokenP,"]")) { *temp++ = *tokenP++; *tokenType = RIGHTCOMMENT; } else if (IsIn(*tokenP,"-")) { *temp++ = *tokenP++; *tokenType = DASH; } else if (IsIn(*tokenP,"?")) { *temp++ = *tokenP++; *tokenType = QUESTIONMARK; } else if (IsIn(*tokenP,"abcdefghijklmnopqrstuvwxyzABCDEFGHIJKLMNOPQRSTUVWXYZ_0123456789.")) { allNumbers = TRUE; if (!IsIn(*tokenP,"0123456789.")) allNumbers = FALSE; *temp++ = *tokenP++; while(IsIn(*tokenP,"abcdefghijklmnopqrstuvwxyzABCDEFGHIJKLMNOPQRSTUVWXYZ_0123456789.")) { if (!IsIn(*tokenP,"0123456789.")) allNumbers = FALSE; *temp++ = *tokenP++; } if (allNumbers == TRUE) *tokenType = NUMBER; else *tokenType = ALPHA; } else if (IsIn(*tokenP,"*")) { *temp++ = *tokenP++; *tokenType = ASTERISK; } else if (IsIn(*tokenP,"/")) { *temp++ = *tokenP++; *tokenType = BACKSLASH; } else if (IsIn(*tokenP,"'\\'")) { *temp++ = *tokenP++; *tokenType = FORWARDSLASH; } else if (IsIn(*tokenP,"!")) { *temp++ = *tokenP++; *tokenType = EXCLAMATIONMARK; } else if (IsIn(*tokenP,"%")) { *temp++ = *tokenP++; *tokenType = PERCENT; } else if (IsIn(*tokenP,"&")) { *temp++ = *tokenP++; *tokenType = AMPERSTAND; } else if (IsIn(*tokenP,"~")) { *temp++ = *tokenP++; *tokenType = TILDE; } else if (IsIn(*tokenP,"+")) { *temp++ = *tokenP++; *tokenType = PLUS; } else if (IsIn(*tokenP,"^")) { *temp++ = *tokenP++; *tokenType = CARAT; } else if (IsIn(*tokenP,"$")) { *temp++ = *tokenP++; *tokenType = DOLLARSIGN; } else if (IsIn(*tokenP,"@")) { *temp++ = *tokenP++; *tokenType = ATSIGN; } else if (IsIn(*tokenP,"|")) { *temp++ = *tokenP++; *tokenType = PIPESIGN; } else if (IsIn(*tokenP,"{")) { *temp++ = *tokenP++; *tokenType = LEFTCURLY; } else if (IsIn(*tokenP,"}")) { *temp++ = *tokenP++; *tokenType = RIGHTCURLY; } *temp = '\0'; } void GetUserDownPass (TreeNode *p, TreeNode **x, int *y) { if (p != NULL) { GetUserDownPass (p->left, x, y); GetUserDownPass (p->right, x, y); x[(*y)] = p; (*y)++; } } void GetUserHelp (void) { if (generalHelp == YES) { printf ("\n Commands that are available include:\n\n"); printf (" execute -- Executes a file\n"); printf (" quit -- Quits the program\n"); printf (" lset -- Sets the parameters of the likelihood\n"); printf (" prset -- Sets the priors for the parameters\n"); printf (" props -- Sets MCMC chain parameters\n"); printf (" mcmc -- Starts Markov chain Monte Carlo\n"); printf (" mcmcp -- Sets the parameters of a chain (without starting analysis)\n"); printf (" exclude -- Excludes sites from the analysis\n"); printf (" include -- Includes sites\n"); printf (" delete -- Deletes taxa from the analysis\n"); printf (" restore -- Restores taxa\n"); printf (" charstat -- Shows status of characters\n"); printf (" charset -- Assigns a group of sites to a set\n"); printf (" partition -- Assigns a character partition\n"); printf (" constraint -- Assigns a constraint on a tree\n"); printf (" calibration -- Assigns a calibration on a tree\n"); printf (" outgroup -- Assigns a single species to be the outgroup\n"); printf (" root -- Roots tree at outgroup species\n"); printf (" deroot -- Deroots tree\n"); printf (" showtree -- Shows the current user tree\n"); printf (" shownodes -- Shows the current user tree as linked nodes\n"); printf (" sumt -- Summarize the trees in a tree file\n"); printf (" set -- Sets run conditions of program\n"); printf (" usertree -- Assigns a user tree\n"); printf ("\n"); } else { if (!strcmp(helpRequest, "execute")) { printf ("\n execute\n\n"); printf (" This command executes a file called .\n"); printf (" The correct usage is: \"execute \"\n"); printf (" For example, \"execute replicase.nex\" would execute the\n"); printf (" file named \"replicase.nex\". This file must be in the\n"); printf (" same directory as the executable.\n\n"); } else if (!strcmp(helpRequest, "quit")) { printf ("\n quit\n\n"); printf (" This command quits the program.\n"); printf (" The correct usage is: \"quit\"\n\n"); } else if (!strcmp(helpRequest, "lset")) { printf ("\n lset\n\n"); printf (" This command sets the parameters of the likelihood model.\n"); printf (" The likelihood function is the probability of observing the\n"); printf (" data conditional on the phylogenetic model. In order to\n"); printf (" calculate the likelihood, you must assume a model of character\n"); printf (" change. This command lets you tailor these assumptions.\n"); printf (" The correct usage is: \"lset =