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