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