tree.cpp

   1 /*
   2  *  tree.cpp
   3  *  Mothur
   4  *
   5  *  Created by Sarah Westcott on 1/22/09.
   6  *  Copyright 2009 Schloss Lab UMASS Amherst. All rights reserved.
   7  *
   8  */
   9
  10 #include "tree.h"
  11
  12 /*****************************************************************/
  13 Tree::Tree(string g) {
  14         try {
  15                 globaldata = GlobalData::getInstance();
  16                 m = MothurOut::getInstance();
  17
  18                 parseTreeFile();  globaldata->runParse = false;
  19         }
  20         catch(exception& e) {
  21                 m->errorOut(e, "Tree", "Tree - just parse");
  22                 exit(1);
  23         }
  24 }
  25 /*****************************************************************/
  26 Tree::Tree() {
  27         try {
  28                 globaldata = GlobalData::getInstance();
  29                 m = MothurOut::getInstance();
  30
  31                 if (globaldata->runParse == true) {  parseTreeFile();  globaldata->runParse = false;  }
  32 //for(int i = 0; i <    globaldata->Treenames.size(); i++) { cout << i << '\t' << globaldata->Treenames[i] << endl;  }
  33                 numLeaves = globaldata->Treenames.size();
  34                 numNodes = 2*numLeaves - 1;
  35
  36                 tree.resize(numNodes);
  37
  38                 //initialize tree with correct number of nodes, name and group info.
  39                 for (int i = 0; i < numNodes; i++) {
  40                         //initialize leaf nodes
  41                         if (i <= (numLeaves-1)) {
  42                                 tree[i].setName(globaldata->Treenames[i]);
  43                                 vector<string> tempGroups; tempGroups.push_back(globaldata->gTreemap->getGroup(globaldata->Treenames[i]));
  44                                 tree[i].setGroup(tempGroups);
  45                                 //set pcount and pGroup for groupname to 1.
  46                                 tree[i].pcount[globaldata->gTreemap->getGroup(globaldata->Treenames[i])] = 1;
  47                                 tree[i].pGroups[globaldata->gTreemap->getGroup(globaldata->Treenames[i])] = 1;
  48                                 //Treemap knows name, group and index to speed up search
  49                                 globaldata->gTreemap->setIndex(globaldata->Treenames[i], i);
  50
  51                         //intialize non leaf nodes
  52                         }else if (i > (numLeaves-1)) {
  53                                 tree[i].setName("");
  54                                 vector<string> tempGroups;
  55                                 tree[i].setGroup(tempGroups);
  56                         }
  57                 }
  58         }
  59         catch(exception& e) {
  60                 m->errorOut(e, "Tree", "Tree");
  61                 exit(1);
  62         }
  63 }
  64
  65 /*****************************************************************/
  66 Tree::~Tree() {}
  67 /*****************************************************************/
  68 void Tree::addNamesToCounts() {
  69         try {
  70                 //ex. seq1      seq2,seq3,se4
  71                 //              seq1 = pasture
  72                 //              seq2 = forest
  73                 //              seq4 = pasture
  74                 //              seq3 = ocean
  75
  76                 //before this function seq1.pcount = pasture -> 1
  77                 //after                            seq1.pcount = pasture -> 2, forest -> 1, ocean -> 1
  78
  79                 //before this function seq1.pgroups = pasture -> 1
  80                 //after                            seq1.pgroups = pasture -> 1 since that is the dominant group
  81
  82
  83                 //go through each leaf and update its pcounts and pgroups
  84
  85                 //float A = clock();
  86
  87                 for (int i = 0; i < numLeaves; i++) {
  88
  89                         string name = tree[i].getName();
  90
  91                         map<string, string>::iterator itNames = globaldata->names.find(name);
  92
  93                         if (itNames == globaldata->names.end()) { m->mothurOut(name + " is not in your name file, please correct."); m->mothurOutEndLine(); exit(1);  }
  94                         else {
  95                                 vector<string> dupNames;
  96                                 m->splitAtComma(globaldata->names[name], dupNames);
  97
  98                                 map<string, int>::iterator itCounts;
  99                                 int maxPars = 1;
 100                                 for (int j = 0; j < dupNames.size(); j++) {
 101
 102                                         if (dupNames[j] != name) {//you already added yourself in the constructor
 103                                                 string group = globaldata->gTreemap->getGroup(dupNames[j]);
 104
 105                                                 //update pcounts
 106                                                 itCounts = tree[i].pcount.find(group);
 107                                                 if (itCounts == tree[i].pcount.end()) { //new group, add it
 108                                                         tree[i].pcount[group] = 1;
 109                                                 }else {
 110                                                         tree[i].pcount[group]++;
 111                                                 }
 112
 113                                                 //update pgroups
 114                                                 itCounts = tree[i].pGroups.find(group);
 115                                                 if (itCounts == tree[i].pGroups.end()) { //new group, add it
 116                                                         tree[i].pGroups[group] = 1;
 117                                                 }else {
 118                                                         tree[i].pGroups[group]++;
 119                                                 }
 120
 121                                                 //keep highest group
 122                                                 if(tree[i].pGroups[group] > maxPars){
 123                                                         maxPars = tree[i].pGroups[group];
 124                                                 }
 125                                         }//end if
 126                                 }//end for
 127
 128                                 if (maxPars > 1) { //then we have some more dominant groups
 129                                         //erase all the groups that are less than maxPars because you found a more dominant group.
 130                                         for(it=tree[i].pGroups.begin();it!=tree[i].pGroups.end();){
 131                                                 if(it->second < maxPars){
 132                                                         tree[i].pGroups.erase(it++);
 133                                                 }else { it++; }
 134                                         }
 135                                         //set one remaining groups to 1
 136                                         for(it=tree[i].pGroups.begin();it!=tree[i].pGroups.end();it++){
 137                                                 tree[i].pGroups[it->first] = 1;
 138                                         }
 139                                 }//end if
 140
 141                                 //update groups to reflect all the groups this node represents
 142                                 vector<string> nodeGroups;
 143                                 map<string, int>::iterator itGroups;
 144                                 for (itGroups = tree[i].pcount.begin(); itGroups != tree[i].pcount.end(); itGroups++) {
 145                                         nodeGroups.push_back(itGroups->first);
 146                                 }
 147                                 tree[i].setGroup(nodeGroups);
 148
 149                         }//end else
 150                 }//end for
 151
 152                 //float B = clock();
 153                 //cout << "addNamesToCounts\t" << (B - A) / CLOCKS_PER_SEC << endl;
 154
 155         }
 156         catch(exception& e) {
 157                 m->errorOut(e, "Tree", "addNamesToCounts");
 158                 exit(1);
 159         }
 160 }
 161 /*****************************************************************/
 162 int Tree::getIndex(string searchName) {
 163         try {
 164                 //Treemap knows name, group and index to speed up search
 165                 // getIndex function will return the vector index or -1 if seq is not found.
 166                 int index = globaldata->gTreemap->getIndex(searchName);
 167                 return index;
 168
 169         }
 170         catch(exception& e) {
 171                 m->errorOut(e, "Tree", "getIndex");
 172                 exit(1);
 173         }
 174 }
 175 /*****************************************************************/
 176
 177 void Tree::setIndex(string searchName, int index) {
 178         try {
 179                 //set index in treemap
 180                 globaldata->gTreemap->setIndex(searchName, index);
 181         }
 182         catch(exception& e) {
 183                 m->errorOut(e, "Tree", "setIndex");
 184                 exit(1);
 185         }
 186 }
 187 /*****************************************************************/
 188 int Tree::assembleTree() {
 189         try {
 190                 //float A = clock();
 191
 192                 //if user has given a names file we want to include that info in the pgroups and pcount info.
 193                 if(globaldata->names.size() != 0) {  addNamesToCounts();  }
 194
 195                 //build the pGroups in non leaf nodes to be used in the parsimony calcs.
 196                 for (int i = numLeaves; i < numNodes; i++) {
 197                         if (m->control_pressed) { return 1; }
 198
 199                         tree[i].pGroups = (mergeGroups(i));
 200                         tree[i].pcount = (mergeGcounts(i));
 201                 }
 202                 //float B = clock();
 203                 //cout << "assembleTree\t" << (B-A) / CLOCKS_PER_SEC << endl;
 204                 return 0;
 205         }
 206         catch(exception& e) {
 207                 m->errorOut(e, "Tree", "assembleTree");
 208                 exit(1);
 209         }
 210 }
 211 /*****************************************************************/
 212 void Tree::getCopy(Tree* copy) {
 213         try {
 214
 215                 //for each node in the tree copy its info
 216                 for (int i = 0; i < numNodes; i++) {
 217                         //copy name
 218                         tree[i].setName(copy->tree[i].getName());
 219
 220                         //copy group
 221                         tree[i].setGroup(copy->tree[i].getGroup());
 222
 223                         //copy branch length
 224                         tree[i].setBranchLength(copy->tree[i].getBranchLength());
 225
 226                         //copy parent
 227                         tree[i].setParent(copy->tree[i].getParent());
 228
 229                         //copy children
 230                         tree[i].setChildren(copy->tree[i].getLChild(), copy->tree[i].getRChild());
 231
 232                         //copy index in node and tmap
 233                         tree[i].setIndex(copy->tree[i].getIndex());
 234                         setIndex(copy->tree[i].getName(), getIndex(copy->tree[i].getName()));
 235
 236                         //copy pGroups
 237                         tree[i].pGroups = copy->tree[i].pGroups;
 238
 239                         //copy pcount
 240                         tree[i].pcount = copy->tree[i].pcount;
 241                 }
 242
 243         }
 244         catch(exception& e) {
 245                 m->errorOut(e, "Tree", "getCopy");
 246                 exit(1);
 247         }
 248 }
 249 /*****************************************************************/
 250 //returns a map with a groupname and the number of times that group was seen in the children
 251 //for instance if your children are white and black then it would return a map with 2 entries
 252 // p[white] = 1 and p[black] = 1.  Now go up a level and merge that with a node who has p[white] = 1
 253 //and you get p[white] = 2, p[black] = 1, but you erase the p[black] because you have a p value higher than 1.
 254
 255 map<string, int> Tree::mergeGroups(int i) {
 256         try {
 257                 int lc = tree[i].getLChild();
 258                 int rc = tree[i].getRChild();
 259
 260                 //set parsimony groups to left child
 261                 map<string,int> parsimony = tree[lc].pGroups;
 262
 263                 int maxPars = 1;
 264
 265                 //look at right child groups and update maxPars if right child has something higher for that group.
 266                 for(it=tree[rc].pGroups.begin();it!=tree[rc].pGroups.end();it++){
 267                         it2 = parsimony.find(it->first);
 268                         if (it2 != parsimony.end()) {
 269                                 parsimony[it->first]++;
 270                         }else {
 271                                 parsimony[it->first] = 1;
 272                         }
 273
 274                         if(parsimony[it->first] > maxPars){
 275                                 maxPars = parsimony[it->first];
 276                         }
 277                 }
 278
 279                 // this is true if right child had a greater parsimony for a certain group
 280                 if(maxPars > 1){
 281                         //erase all the groups that are only 1 because you found something with 2.
 282                         for(it=parsimony.begin();it!=parsimony.end();){
 283                                 if(it->second == 1){
 284                                         parsimony.erase(it++);
 285                                 }else { it++; }
 286                         }
 287                         //set one remaining groups to 1
 288                         //so with our above example p[white] = 2 would be left and it would become p[white] = 1
 289                         for(it=parsimony.begin();it!=parsimony.end();it++){
 290                                 parsimony[it->first] = 1;
 291                         }
 292
 293                 }
 294
 295                 return parsimony;
 296         }
 297         catch(exception& e) {
 298                 m->errorOut(e, "Tree", "mergeGroups");
 299                 exit(1);
 300         }
 301 }
 302 /*****************************************************************/
 303 //returns a map with a groupname and the number of times that group was seen in the children
 304 //for instance if your children are white and black then it would return a map with 2 entries
 305 // p[white] = 1 and p[black] = 1.  Now go up a level and merge that with a node who has p[white] = 1
 306 //and you get p[white] = 2, p[black] = 1, but you erase the p[black] because you have a p value higher than 1.
 307
 308 map<string, int> Tree::mergeUserGroups(int i, vector<string> g) {
 309         try {
 310
 311                 int lc = tree[i].getLChild();
 312                 int rc = tree[i].getRChild();
 313
 314                 //loop through nodes groups removing the ones the user doesn't want
 315                 for(it=tree[lc].pGroups.begin();it!=tree[lc].pGroups.end();){
 316                                 if (m->inUsersGroups(it->first, g) != true) {
 317                                         tree[lc].pGroups.erase(it++);
 318                                 }else { it++; }
 319                 }
 320
 321                 //loop through nodes groups removing the ones the user doesn't want
 322                 for(it=tree[rc].pGroups.begin();it!=tree[rc].pGroups.end();){
 323                                 if (m->inUsersGroups(it->first, g) != true) {
 324                                         tree[rc].pGroups.erase(it++);
 325                                 }else { it++; }
 326                 }
 327
 328                 //set parsimony groups to left child
 329                 map<string,int> parsimony = tree[lc].pGroups;
 330
 331                 int maxPars = 1;
 332
 333                 //look at right child groups and update maxPars if right child has something higher for that group.
 334                 for(it=tree[rc].pGroups.begin();it!=tree[rc].pGroups.end();it++){
 335                         it2 = parsimony.find(it->first);
 336                         if (it2 != parsimony.end()) {
 337                                 parsimony[it->first]++;
 338                         }else {
 339                                 parsimony[it->first] = 1;
 340                         }
 341
 342                         if(parsimony[it->first] > maxPars){
 343                                 maxPars = parsimony[it->first];
 344                         }
 345                 }
 346
 347                 // this is true if right child had a greater parsimony for a certain group
 348                 if(maxPars > 1){
 349                         //erase all the groups that are only 1 because you found something with 2.
 350                         for(it=parsimony.begin();it!=parsimony.end();){
 351                                 if(it->second == 1){
 352                                         parsimony.erase(it++);
 353                                 }else { it++; }
 354                         }
 355
 356                         for(it=parsimony.begin();it!=parsimony.end();it++){
 357                                 parsimony[it->first] = 1;
 358                         }
 359                 }
 360
 361                 return parsimony;
 362         }
 363         catch(exception& e) {
 364                 m->errorOut(e, "Tree", "mergeUserGroups");
 365                 exit(1);
 366         }
 367 }
 368
 369
 370 /**************************************************************************************************/
 371
 372 map<string,int> Tree::mergeGcounts(int position) {
 373         try{
 374                 map<string,int>::iterator pos;
 375
 376                 int lc = tree[position].getLChild();
 377                 int rc = tree[position].getRChild();
 378
 379                 map<string,int> sum = tree[lc].pcount;
 380
 381                 for(it=tree[rc].pcount.begin();it!=tree[rc].pcount.end();it++){
 382                         sum[it->first] += it->second;
 383                 }
 384                 return sum;
 385         }
 386         catch(exception& e) {
 387                 m->errorOut(e, "Tree", "mergeGcounts");
 388                 exit(1);
 389         }
 390 }
 391 /**************************************************************************************************/
 392
 393 void Tree::randomLabels(vector<string> g) {
 394         try {
 395
 396                 for(int i = 0; i < numLeaves; i++){
 397                         int z;
 398                         //get random index to switch with
 399                         z = int((float)(i+1) * (float)(rand()) / ((float)RAND_MAX+1.0));
 400
 401                         //you only want to randomize the nodes that are from a group the user wants analyzed, so
 402                         //if either of the leaf nodes you are about to switch are not in the users groups then you don't want to switch them.
 403                         bool treez, treei;
 404
 405                         treez = m->inUsersGroups(tree[z].getGroup(), g);
 406                         treei = m->inUsersGroups(tree[i].getGroup(), g);
 407
 408                         if ((treez == true) && (treei == true)) {
 409                                 //switches node i and node z's info.
 410                                 map<string,int> lib_hold = tree[z].pGroups;
 411                                 tree[z].pGroups = (tree[i].pGroups);
 412                                 tree[i].pGroups = (lib_hold);
 413
 414                                 vector<string> zgroup = tree[z].getGroup();
 415                                 tree[z].setGroup(tree[i].getGroup());
 416                                 tree[i].setGroup(zgroup);
 417
 418                                 string zname = tree[z].getName();
 419                                 tree[z].setName(tree[i].getName());
 420                                 tree[i].setName(zname);
 421
 422                                 map<string,int> gcount_hold = tree[z].pcount;
 423                                 tree[z].pcount = (tree[i].pcount);
 424                                 tree[i].pcount = (gcount_hold);
 425                         }
 426                 }
 427         }
 428         catch(exception& e) {
 429                 m->errorOut(e, "Tree", "randomLabels");
 430                 exit(1);
 431         }
 432 }
 433 /**************************************************************************************************
 434
 435 void Tree::randomLabels(string groupA, string groupB) {
 436         try {
 437                 int numSeqsA = globaldata->gTreemap->seqsPerGroup[groupA];
 438                 int numSeqsB = globaldata->gTreemap->seqsPerGroup[groupB];
 439
 440                 vector<string> randomGroups(numSeqsA+numSeqsB, groupA);
 441                 for(int i=numSeqsA;i<randomGroups.size();i++){
 442                         randomGroups[i] = groupB;
 443                 }
 444                 random_shuffle(randomGroups.begin(), randomGroups.end());
 445
 446                 int randomCounter = 0;
 447                 for(int i=0;i<numLeaves;i++){
 448                         if(tree[i].getGroup() == groupA || tree[i].getGroup() == groupB){
 449                                 tree[i].setGroup(randomGroups[randomCounter]);
 450                                 tree[i].pcount.clear();
 451                                 tree[i].pcount[randomGroups[randomCounter]] = 1;
 452                                 tree[i].pGroups.clear();
 453                                 tree[i].pGroups[randomGroups[randomCounter]] = 1;
 454                                 randomCounter++;
 455                         }
 456                 }
 457         }
 458         catch(exception& e) {
 459                 m->errorOut(e, "Tree", "randomLabels");
 460                 exit(1);
 461         }
 462 }
 463 /**************************************************************************************************/
 464 void Tree::randomBlengths()  {
 465         try {
 466                 for(int i=numNodes-1;i>=0;i--){
 467                         int z = int((float)(i+1) * (float)(rand()) / ((float)RAND_MAX+1.0));
 468
 469                         float bl_hold = tree[z].getBranchLength();
 470                         tree[z].setBranchLength(tree[i].getBranchLength());
 471                         tree[i].setBranchLength(bl_hold);
 472                 }
 473         }
 474         catch(exception& e) {
 475                 m->errorOut(e, "Tree", "randomBlengths");
 476                 exit(1);
 477         }
 478 }
 479 /*************************************************************************************************/
 480 void Tree::assembleRandomUnifracTree(vector<string> g) {
 481         randomLabels(g);
 482         assembleTree();
 483 }
 484 /*************************************************************************************************/
 485 void Tree::assembleRandomUnifracTree(string groupA, string groupB) {
 486
 487         vector<string> temp; temp.push_back(groupA); temp.push_back(groupB);
 488         randomLabels(temp);
 489         assembleTree();
 490 }
 491
 492 /*************************************************************************************************/
 493 //for now it's just random topology but may become random labels as well later that why this is such a simple function now...
 494 void Tree::assembleRandomTree() {
 495         randomTopology();
 496         assembleTree();
 497 }
 498 /**************************************************************************************************/
 499
 500 void Tree::randomTopology() {
 501         try {
 502                 for(int i=0;i<numNodes;i++){
 503                         tree[i].setParent(-1);
 504                 }
 505                 for(int i=numLeaves;i<numNodes;i++){
 506                         tree[i].setChildren(-1, -1);
 507                 }
 508
 509                 for(int i=numLeaves;i<numNodes;i++){
 510                         int escape =0;
 511                         int rnd_index1, rnd_index2;
 512                         while(escape == 0){
 513                                 rnd_index1 = (int)(((double)rand() / (double) RAND_MAX)*i);
 514                                 if(tree[rnd_index1].getParent() == -1){escape = 1;}
 515                         }
 516
 517                         escape = 0;
 518                         while(escape == 0){
 519                                 rnd_index2 = (int)(((double)rand() / (double) RAND_MAX)*i);
 520                                 if(rnd_index2 != rnd_index1 && tree[rnd_index2].getParent() == -1){
 521                                         escape = 1;
 522                                 }
 523                         }
 524
 525                         tree[i].setChildren(rnd_index1,rnd_index2);
 526                         tree[i].setParent(-1);
 527                         tree[rnd_index1].setParent(i);
 528                         tree[rnd_index2].setParent(i);
 529                 }
 530         }
 531         catch(exception& e) {
 532                 m->errorOut(e, "Tree", "randomTopology");
 533                 exit(1);
 534         }
 535 }
 536 /*****************************************************************/
 537 void Tree::print(ostream& out) {
 538         try {
 539                 int root = findRoot();
 540                 printBranch(root, out, "branch");
 541                 out << ";" << endl;
 542         }
 543         catch(exception& e) {
 544                 m->errorOut(e, "Tree", "print");
 545                 exit(1);
 546         }
 547 }
 548 /*****************************************************************/
 549 void Tree::printForBoot(ostream& out) {
 550         try {
 551                 int root = findRoot();
 552                 printBranch(root, out, "boot");
 553                 out << ";" << endl;
 554         }
 555         catch(exception& e) {
 556                 m->errorOut(e, "Tree", "printForBoot");
 557                 exit(1);
 558         }
 559 }
 560
 561 /*****************************************************************/
 562 // This prints out the tree in Newick form.
 563 void Tree::createNewickFile(string f) {
 564         try {
 565                 int root = findRoot();
 566                 //filename = m->getRootName(globaldata->getTreeFile()) + "newick";
 567                 filename = f;
 568
 569                 m->openOutputFile(filename, out);
 570
 571                 printBranch(root, out, "branch");
 572
 573                 // you are at the end of the tree
 574                 out << ";" << endl;
 575                 out.close();
 576         }
 577         catch(exception& e) {
 578                 m->errorOut(e, "Tree", "createNewickFile");
 579                 exit(1);
 580         }
 581 }
 582
 583 /*****************************************************************/
 584 //This function finds the index of the root node.
 585
 586 int Tree::findRoot() {
 587         try {
 588                 for (int i = 0; i < numNodes; i++) {
 589                         //you found the root
 590                         if (tree[i].getParent() == -1) { return i; }
 591                         //cout << "i = " << i << endl;
 592                         //cout << "i's parent = " << tree[i].getParent() << endl;
 593                 }
 594                 return -1;
 595         }
 596         catch(exception& e) {
 597                 m->errorOut(e, "Tree", "findRoot");
 598                 exit(1);
 599         }
 600 }
 601
 602 /*****************************************************************/
 603 void Tree::printBranch(int node, ostream& out, string mode) {
 604         try {
 605
 606                 // you are not a leaf
 607                 if (tree[node].getLChild() != -1) {
 608                         out << "(";
 609                         printBranch(tree[node].getLChild(), out, mode);
 610                         out << ",";
 611                         printBranch(tree[node].getRChild(), out, mode);
 612                         out << ")";
 613                         if (mode == "branch") {
 614                                 //if there is a branch length then print it
 615                                 if (tree[node].getBranchLength() != -1) {
 616                                         out << ":" << tree[node].getBranchLength();
 617                                 }
 618                         }else if (mode == "boot") {
 619                                 //if there is a label then print it
 620                                 if (tree[node].getLabel() != -1) {
 621                                         out << tree[node].getLabel();
 622                                 }
 623                         }
 624                 }else { //you are a leaf
 625                         string leafGroup = globaldata->gTreemap->getGroup(tree[node].getName());
 626
 627                         out << leafGroup;
 628                         if (mode == "branch") {
 629                                 //if there is a branch length then print it
 630                                 if (tree[node].getBranchLength() != -1) {
 631                                         out << ":" << tree[node].getBranchLength();
 632                                 }
 633                         }else if (mode == "boot") {
 634                                 //if there is a label then print it
 635                                 if (tree[node].getLabel() != -1) {
 636                                         out << tree[node].getLabel();
 637                                 }
 638                         }
 639                 }
 640
 641         }
 642         catch(exception& e) {
 643                 m->errorOut(e, "Tree", "printBranch");
 644                 exit(1);
 645         }
 646 }
 647
 648 /*****************************************************************/
 649
 650 void Tree::printTree() {
 651
 652         for(int i=0;i<numNodes;i++){
 653                 cout << i << '\t';
 654                 tree[i].printNode();
 655         }
 656
 657 }
 658
 659 /*****************************************************************/
 660 //this code is a mess and should be rethought...-slw
 661 void Tree::parseTreeFile() {
 662
 663         //only takes names from the first tree and assumes that all trees use the same names.
 664         try {
 665                 string filename = globaldata->getTreeFile();
 666                 ifstream filehandle;
 667                 m->openInputFile(filename, filehandle);
 668                 int c, comment;
 669                 comment = 0;
 670                 int done = 1;
 671
 672                 //ifyou are not a nexus file
 673                 if((c = filehandle.peek()) != '#') {
 674                         while((c = filehandle.peek()) != ';') {
 675                                 while ((c = filehandle.peek()) != ';') {
 676                                         // get past comments
 677                                         if(c == '[') {
 678                                                 comment = 1;
 679                                         }
 680                                         if(c == ']'){
 681                                                 comment = 0;
 682                                         }
 683                                         if((c == '(') && (comment != 1)){ break; }
 684                                         filehandle.get();
 685                                 }
 686
 687                                 done = readTreeString(filehandle);
 688                                 if (done == 0) { break; }
 689                         }
 690                 //ifyou are a nexus file
 691                 }else if((c = filehandle.peek()) == '#') {
 692                         string holder = "";
 693
 694                         // get past comments
 695                         while(holder != "translate" && holder != "Translate"){
 696                                 if(holder == "[" || holder == "[!"){
 697                                         comment = 1;
 698                                 }
 699                                 if(holder == "]"){
 700                                         comment = 0;
 701                                 }
 702                                 filehandle >> holder;
 703
 704                                 //if there is no translate then you must read tree string otherwise use translate to get names
 705                                 if((holder == "tree") && (comment != 1)){
 706                                         //pass over the "tree rep.6878900 = "
 707                                         while (((c = filehandle.get()) != '(') && ((c = filehandle.peek()) != EOF)) {;}
 708
 709                                         if(c == EOF) { break; }
 710                                         filehandle.putback(c);  //put back first ( of tree.
 711                                         done = readTreeString(filehandle);
 712
 713                                         break;
 714                                 }
 715
 716                                 if (done == 0) { break;  }
 717                         }
 718
 719                         //use nexus translation rather than parsing tree to save time
 720                         if((holder == "translate") || (holder == "Translate")) {
 721
 722                                 string number, name, h;
 723                                 h = ""; // so it enters the loop the first time
 724                                 while((h != ";") && (number != ";")) {
 725                                         filehandle >> number;
 726                                         filehandle >> name;
 727
 728                                         //c = , until done with translation then c = ;
 729                                         h = name.substr(name.length()-1, name.length());
 730                                         name.erase(name.end()-1);  //erase the comma
 731                                         globaldata->Treenames.push_back(number);
 732                                 }
 733                                 if(number == ";") { globaldata->Treenames.pop_back(); }  //in case ';' from translation is on next line instead of next to last name
 734                         }
 735                 }
 736                 filehandle.close();
 737
 738                 //for (int i = 0; i < globaldata->Treenames.size(); i++) {
 739 //cout << globaldata->Treenames[i] << endl; }
 740 //cout << globaldata->Treenames.size() << endl;
 741         }
 742         catch(exception& e) {
 743                 m->errorOut(e, "Tree", "parseTreeFile");
 744                 exit(1);
 745         }
 746 }
 747 /*******************************************************/
 748
 749 /*******************************************************/
 750 int Tree::readTreeString(ifstream& filehandle)  {
 751         try {
 752                 int c;
 753                 string name;  //, k
 754
 755                 while((c = filehandle.peek()) != ';') {
 756 //k = c;
 757 //cout << " at beginning of while " <<  k << endl;
 758                         if(c == ')')  {
 759                                 //to pass over labels in trees
 760                                 c=filehandle.get();
 761                                 while((c!=',') && (c != -1) && (c!= ':') && (c!=';')){ c=filehandle.get(); }
 762                                 filehandle.putback(c);
 763                         }
 764                         if(c == ';') { return 0; }
 765                         if(c == -1) { return 0; }
 766                         //if you are a name
 767                         if((c != '(') && (c != ')') && (c != ',') && (c != ':') && (c != '\n') && (c != '\t') && (c != 32)) { //32 is space
 768                                 name = "";
 769                                 c = filehandle.get();
 770                         //k = c;
 771 //cout << k << endl;
 772                                 while ((c != '(') && (c != ')') && (c != ',') && (c != ':') && (c != '\n') && (c != 32) && (c != '\t')) {
 773                                         name += c;
 774                                         c = filehandle.get();
 775                         //k = c;
 776 //cout << " in name while " << k << endl;
 777                                 }
 778
 779 //cout << "name = " << name << endl;
 780                                 globaldata->Treenames.push_back(name);
 781                                 filehandle.putback(c);
 782 //k = c;
 783 //cout << " after putback" <<  k << endl;
 784                         }
 785
 786                         if(c  == ':') { //read until you reach the end of the branch length
 787                                 while ((c != '(') && (c != ')') && (c != ',') && (c != ';') && (c != '\n') && (c != '\t') && (c != 32)) {
 788                                         c = filehandle.get();
 789         //k = c;
 790         //cout << " in branch while " << k << endl;
 791                                 }
 792                                 filehandle.putback(c);
 793                         }
 794
 795                         c = filehandle.get();
 796 //k = c;
 797         //cout << " here after get " << k << endl;
 798                         if(c == ';') { return 0; }
 799                         if(c == ')') { filehandle.putback(c); }
 800         //k = c;
 801 //cout << k << endl;
 802
 803                 }
 804                 return 0;
 805         }
 806         catch(exception& e) {
 807                 m->errorOut(e, "Tree", "readTreeString");
 808                 exit(1);
 809         }
 810 }
 811
 812 /*******************************************************/
 813
 814 /*******************************************************/
 815