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(int num, CountTable* t) : ct(t) {
  14         try {
  15                 m = MothurOut::getInstance();
  16
  17                 numLeaves = num;
  18                 numNodes = 2*numLeaves - 1;
  19
  20                 tree.resize(numNodes);
  21         }
  22         catch(exception& e) {
  23                 m->errorOut(e, "Tree", "Tree - numNodes");
  24                 exit(1);
  25         }
  26 }
  27 /*****************************************************************/
  28 Tree::Tree(string g) { //do not use tree generated by this its just to extract the treenames, its a chicken before the egg thing that needs to be revisited.
  29         try {
  30                 m = MothurOut::getInstance();
  31                 parseTreeFile();  m->runParse = false;
  32         }
  33         catch(exception& e) {
  34                 m->errorOut(e, "Tree", "Tree - just parse");
  35                 exit(1);
  36         }
  37 }
  38 /*****************************************************************/
  39 Tree::Tree(CountTable* t) : ct(t) {
  40         try {
  41                 m = MothurOut::getInstance();
  42
  43                 if (m->runParse == true) {  parseTreeFile();  m->runParse = false;  }
  44
  45                 numLeaves = m->Treenames.size();
  46                 numNodes = 2*numLeaves - 1;
  47
  48                 tree.resize(numNodes);
  49
  50                 //initialize groupNodeInfo
  51         vector<string> namesOfGroups = ct->getNamesOfGroups();
  52                 for (int i = 0; i < namesOfGroups.size(); i++) {  groupNodeInfo[namesOfGroups[i]].resize(0);  }
  53
  54                 //initialize tree with correct number of nodes, name and group info.
  55                 for (int i = 0; i < numNodes; i++) {
  56                         //initialize leaf nodes
  57                         if (i <= (numLeaves-1)) {
  58                                 tree[i].setName(m->Treenames[i]);
  59
  60                                 //save group info
  61                 int maxPars = 1;
  62                                 vector<string> group;
  63                 vector<int> counts = ct->getGroupCounts(m->Treenames[i]);
  64                                 for (int j = 0; j < namesOfGroups.size(); j++) {
  65                     if (counts[j] != 0) { //you have seqs from this group
  66                         groupNodeInfo[namesOfGroups[j]].push_back(i);
  67                         group.push_back(namesOfGroups[j]);
  68                         tree[i].pGroups[namesOfGroups[j]] = counts[j];
  69                         tree[i].pcount[namesOfGroups[j]] = counts[j];
  70                         //keep highest group
  71                                                 if(counts[j] > maxPars){ maxPars = counts[j]; }
  72                     }
  73                 }
  74                                 tree[i].setGroup(group);
  75                                 setIndex(m->Treenames[i], i);
  76
  77                 if (maxPars > 1) { //then we have some more dominant groups
  78                                         //erase all the groups that are less than maxPars because you found a more dominant group.
  79                                         for(it=tree[i].pGroups.begin();it!=tree[i].pGroups.end();){
  80                                                 if(it->second < maxPars){
  81                                                         tree[i].pGroups.erase(it++);
  82                                                 }else { it++; }
  83                                         }
  84                                         //set one remaining groups to 1
  85                                         for(it=tree[i].pGroups.begin();it!=tree[i].pGroups.end();it++){
  86                                                 tree[i].pGroups[it->first] = 1;
  87                                         }
  88                                 }//end if
  89
  90                         //intialize non leaf nodes
  91                         }else if (i > (numLeaves-1)) {
  92                                 tree[i].setName("");
  93                                 vector<string> tempGroups;
  94                                 tree[i].setGroup(tempGroups);
  95                         }
  96                 }
  97
  98         }
  99         catch(exception& e) {
 100                 m->errorOut(e, "Tree", "Tree");
 101                 exit(1);
 102         }
 103 }
 104 /*****************************************************************/
 105 Tree::Tree(CountTable* t, vector< vector<double> >& sims) : ct(t) {
 106         try {
 107                 m = MothurOut::getInstance();
 108
 109                 if (m->runParse == true) {  parseTreeFile();  m->runParse = false;  }
 110                 numLeaves = m->Treenames.size();
 111                 numNodes = 2*numLeaves - 1;
 112
 113                 tree.resize(numNodes);
 114
 115                 //initialize groupNodeInfo
 116         vector<string> namesOfGroups = ct->getNamesOfGroups();
 117                 for (int i = 0; i < namesOfGroups.size(); i++) {  groupNodeInfo[namesOfGroups[i]].resize(0);  }
 118
 119                 //initialize tree with correct number of nodes, name and group info.
 120                 for (int i = 0; i < numNodes; i++) {
 121                         //initialize leaf nodes
 122                         if (i <= (numLeaves-1)) {
 123                                 tree[i].setName(m->Treenames[i]);
 124
 125                                 //save group info
 126                 int maxPars = 1;
 127                                 vector<string> group;
 128                 vector<int> counts = ct->getGroupCounts(m->Treenames[i]);
 129                                 for (int j = 0; j < namesOfGroups.size(); j++) {
 130                     if (counts[j] != 0) { //you have seqs from this group
 131                         groupNodeInfo[namesOfGroups[j]].push_back(i);
 132                         group.push_back(namesOfGroups[j]);
 133                         tree[i].pGroups[namesOfGroups[j]] = counts[j];
 134                         tree[i].pcount[namesOfGroups[j]] = counts[j];
 135                         //keep highest group
 136                                                 if(counts[j] > maxPars){ maxPars = counts[j]; }
 137                     }
 138                 }
 139                                 tree[i].setGroup(group);
 140                                 setIndex(m->Treenames[i], i);
 141
 142                 if (maxPars > 1) { //then we have some more dominant groups
 143                                         //erase all the groups that are less than maxPars because you found a more dominant group.
 144                                         for(it=tree[i].pGroups.begin();it!=tree[i].pGroups.end();){
 145                                                 if(it->second < maxPars){
 146                                                         tree[i].pGroups.erase(it++);
 147                                                 }else { it++; }
 148                                         }
 149                                         //set one remaining groups to 1
 150                                         for(it=tree[i].pGroups.begin();it!=tree[i].pGroups.end();it++){
 151                                                 tree[i].pGroups[it->first] = 1;
 152                                         }
 153                                 }//end if
 154
 155                 //intialize non leaf nodes
 156                         }else if (i > (numLeaves-1)) {
 157                                 tree[i].setName("");
 158                                 vector<string> tempGroups;
 159                                 tree[i].setGroup(tempGroups);
 160                         }
 161                 }
 162
 163
 164         //build tree from matrix
 165         //initialize indexes
 166         map<int, int> thisIndexes;  //maps row in simMatrix to vector index in the tree
 167         for (int g = 0; g < numLeaves; g++) {   thisIndexes[g] = g;     }
 168
 169                 //do merges and create tree structure by setting parents and children
 170                 //there are numGroups - 1 merges to do
 171                 for (int i = 0; i < (numLeaves - 1); i++) {
 172                         float largest = -1000.0;
 173
 174                         if (m->control_pressed) { break; }
 175
 176                         int row, column;
 177                         //find largest value in sims matrix by searching lower triangle
 178                         for (int j = 1; j < sims.size(); j++) {
 179                                 for (int k = 0; k < j; k++) {
 180                                         if (sims[j][k] > largest) {  largest = sims[j][k]; row = j; column = k;  }
 181                                 }
 182                         }
 183
 184                         //set non-leaf node info and update leaves to know their parents
 185                         //non-leaf
 186                         tree[numLeaves + i].setChildren(thisIndexes[row], thisIndexes[column]);
 187
 188                         //parents
 189                         tree[thisIndexes[row]].setParent(numLeaves + i);
 190                         tree[thisIndexes[column]].setParent(numLeaves + i);
 191
 192                         //blength = distance / 2;
 193                         float blength = ((1.0 - largest) / 2);
 194
 195                         //branchlengths
 196                         tree[thisIndexes[row]].setBranchLength(blength - tree[thisIndexes[row]].getLengthToLeaves());
 197                         tree[thisIndexes[column]].setBranchLength(blength - tree[thisIndexes[column]].getLengthToLeaves());
 198
 199                         //set your length to leaves to your childs length plus branchlength
 200                         tree[numLeaves + i].setLengthToLeaves(tree[thisIndexes[row]].getLengthToLeaves() + tree[thisIndexes[row]].getBranchLength());
 201
 202
 203                         //update index
 204                         thisIndexes[row] = numLeaves+i;
 205                         thisIndexes[column] = numLeaves+i;
 206
 207                         //remove highest value that caused the merge.
 208                         sims[row][column] = -1000.0;
 209                         sims[column][row] = -1000.0;
 210
 211                         //merge values in simsMatrix
 212                         for (int n = 0; n < sims.size(); n++)   {
 213                                 //row becomes merge of 2 groups
 214                                 sims[row][n] = (sims[row][n] + sims[column][n]) / 2;
 215                                 sims[n][row] = sims[row][n];
 216                                 //delete column
 217                                 sims[column][n] = -1000.0;
 218                                 sims[n][column] = -1000.0;
 219                         }
 220                 }
 221
 222                 //adjust tree to make sure root to tip length is .5
 223                 int root = findRoot();
 224                 tree[root].setBranchLength((0.5 - tree[root].getLengthToLeaves()));
 225
 226     }
 227         catch(exception& e) {
 228                 m->errorOut(e, "Tree", "Tree");
 229                 exit(1);
 230         }
 231 }
 232 /*****************************************************************/
 233 Tree::~Tree() {}
 234 /*****************************************************************
 235 void Tree::addNamesToCounts(map<string, string> nameMap) {
 236         try {
 237                 //ex. seq1      seq2,seq3,se4
 238                 //              seq1 = pasture
 239                 //              seq2 = forest
 240                 //              seq4 = pasture
 241                 //              seq3 = ocean
 242
 243                 //before this function seq1.pcount = pasture -> 1
 244                 //after                            seq1.pcount = pasture -> 2, forest -> 1, ocean -> 1
 245
 246                 //before this function seq1.pgroups = pasture -> 1
 247                 //after                            seq1.pgroups = pasture -> 1 since that is the dominant group
 248
 249
 250                 //go through each leaf and update its pcounts and pgroups
 251
 252                 //float A = clock();
 253
 254                 for (int i = 0; i < numLeaves; i++) {
 255
 256                         string name = tree[i].getName();
 257
 258                         map<string, string>::iterator itNames = nameMap.find(name);
 259
 260                         if (itNames == nameMap.end()) { m->mothurOut(name + " is not in your name file, please correct."); m->mothurOutEndLine(); exit(1);  }
 261                         else {
 262                                 vector<string> dupNames;
 263                                 m->splitAtComma(nameMap[name], dupNames);
 264
 265                                 map<string, int>::iterator itCounts;
 266                                 int maxPars = 1;
 267                                 set<string> groupsAddedForThisNode;
 268                                 for (int j = 0; j < dupNames.size(); j++) {
 269
 270                                         string group = tmap->getGroup(dupNames[j]);
 271
 272                                         if (dupNames[j] != name) {//you already added yourself in the constructor
 273
 274                                                 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
 275
 276                                                 //update pcounts
 277                                                 itCounts = tree[i].pcount.find(group);
 278                                                 if (itCounts == tree[i].pcount.end()) { //new group, add it
 279                                                         tree[i].pcount[group] = 1;
 280                                                 }else {
 281                                                         tree[i].pcount[group]++;
 282                                                 }
 283
 284                                                 //update pgroups
 285                                                 itCounts = tree[i].pGroups.find(group);
 286                                                 if (itCounts == tree[i].pGroups.end()) { //new group, add it
 287                                                         tree[i].pGroups[group] = 1;
 288                                                 }else{
 289                                                         tree[i].pGroups[group]++;
 290                                                 }
 291
 292                                                 //keep highest group
 293                                                 if(tree[i].pGroups[group] > maxPars){
 294                                                         maxPars = tree[i].pGroups[group];
 295                                                 }
 296                                         }else {  groupsAddedForThisNode.insert(group);  } //add it so you don't add it to groupNodeInfo again
 297                                 }//end for
 298
 299                                 if (maxPars > 1) { //then we have some more dominant groups
 300                                         //erase all the groups that are less than maxPars because you found a more dominant group.
 301                                         for(it=tree[i].pGroups.begin();it!=tree[i].pGroups.end();){
 302                                                 if(it->second < maxPars){
 303                                                         tree[i].pGroups.erase(it++);
 304                                                 }else { it++; }
 305                                         }
 306                                         //set one remaining groups to 1
 307                                         for(it=tree[i].pGroups.begin();it!=tree[i].pGroups.end();it++){
 308                                                 tree[i].pGroups[it->first] = 1;
 309                                         }
 310                                 }//end if
 311
 312                                 //update groups to reflect all the groups this node represents
 313                                 vector<string> nodeGroups;
 314                                 map<string, int>::iterator itGroups;
 315                                 for (itGroups = tree[i].pcount.begin(); itGroups != tree[i].pcount.end(); itGroups++) {
 316                                         nodeGroups.push_back(itGroups->first);
 317                                 }
 318                                 tree[i].setGroup(nodeGroups);
 319
 320                         }//end else
 321                 }//end for
 322
 323                 //float B = clock();
 324                 //cout << "addNamesToCounts\t" << (B - A) / CLOCKS_PER_SEC << endl;
 325
 326         }
 327         catch(exception& e) {
 328                 m->errorOut(e, "Tree", "addNamesToCounts");
 329                 exit(1);
 330         }
 331 }*/
 332 /*****************************************************************/
 333 int Tree::getIndex(string searchName) {
 334         try {
 335         map<string, int>::iterator itIndex = indexes.find(searchName);
 336         if (itIndex != indexes.end()) {
 337             return itIndex->second;
 338         }
 339                 return -1;
 340         }
 341         catch(exception& e) {
 342                 m->errorOut(e, "Tree", "getIndex");
 343                 exit(1);
 344         }
 345 }
 346 /*****************************************************************/
 347
 348 void Tree::setIndex(string searchName, int index) {
 349         try {
 350                 map<string, int>::iterator itIndex = indexes.find(searchName);
 351         if (itIndex == indexes.end()) {
 352             indexes[searchName] = index;
 353         }
 354         }
 355         catch(exception& e) {
 356                 m->errorOut(e, "Tree", "setIndex");
 357                 exit(1);
 358         }
 359 }
 360 /*****************************************************************/
 361 int Tree::assembleTree() {
 362         try {
 363                 //build the pGroups in non leaf nodes to be used in the parsimony calcs.
 364                 for (int i = numLeaves; i < numNodes; i++) {
 365                         if (m->control_pressed) { return 1; }
 366
 367                         tree[i].pGroups = (mergeGroups(i));
 368                         tree[i].pcount = (mergeGcounts(i));
 369                 }
 370
 371                 return 0;
 372         }
 373         catch(exception& e) {
 374                 m->errorOut(e, "Tree", "assembleTree");
 375                 exit(1);
 376         }
 377 }
 378 /*****************************************************************/
 379 //assumes leaf node names are in groups and no names file - used by indicator command
 380 void Tree::getSubTree(Tree* Ctree, vector<string> Groups) {
 381         try {
 382
 383         //copy Tree since we are going to destroy it
 384         Tree* copy = new Tree(ct);
 385         copy->getCopy(Ctree);
 386         copy->assembleTree();
 387
 388                 //we want to select some of the leaf nodes to create the output tree
 389                 //go through the input Tree starting at parents of leaves
 390         //initialize groupNodeInfo
 391         vector<string> namesOfGroups = ct->getNamesOfGroups();
 392                 for (int i = 0; i < namesOfGroups.size(); i++) {  groupNodeInfo[namesOfGroups[i]].resize(0);  }
 393
 394                 //initialize tree with correct number of nodes, name and group info.
 395                 for (int i = 0; i < numNodes; i++) {
 396                         //initialize leaf nodes
 397                         if (i <= (numLeaves-1)) {
 398                                 tree[i].setName(Groups[i]);
 399
 400                                 //save group info
 401                 int maxPars = 1;
 402                                 vector<string> group;
 403                 vector<int> counts = ct->getGroupCounts(Groups[i]);
 404                                 for (int j = 0; j < namesOfGroups.size(); j++) {
 405                     if (counts[j] != 0) { //you have seqs from this group
 406                         groupNodeInfo[namesOfGroups[j]].push_back(i);
 407                         group.push_back(namesOfGroups[j]);
 408                         tree[i].pGroups[namesOfGroups[j]] = counts[j];
 409                         tree[i].pcount[namesOfGroups[j]] = counts[j];
 410                         //keep highest group
 411                                                 if(counts[j] > maxPars){ maxPars = counts[j]; }
 412                     }
 413                 }
 414                                 tree[i].setGroup(group);
 415                                 setIndex(Groups[i], i);
 416
 417                 if (maxPars > 1) { //then we have some more dominant groups
 418                                         //erase all the groups that are less than maxPars because you found a more dominant group.
 419                                         for(it=tree[i].pGroups.begin();it!=tree[i].pGroups.end();){
 420                                                 if(it->second < maxPars){
 421                                                         tree[i].pGroups.erase(it++);
 422                                                 }else { it++; }
 423                                         }
 424                                         //set one remaining groups to 1
 425                                         for(it=tree[i].pGroups.begin();it!=tree[i].pGroups.end();it++){
 426                                                 tree[i].pGroups[it->first] = 1;
 427                                         }
 428                                 }//end if
 429
 430                 //intialize non leaf nodes
 431                         }else if (i > (numLeaves-1)) {
 432                                 tree[i].setName("");
 433                                 vector<string> tempGroups;
 434                                 tree[i].setGroup(tempGroups);
 435                         }
 436                 }
 437
 438                 set<int> removedLeaves;
 439                 for (int i = 0; i < copy->getNumLeaves(); i++) {
 440
 441                         if (removedLeaves.count(i) == 0) {
 442
 443                         //am I in the group
 444                         int parent = copy->tree[i].getParent();
 445
 446                         if (parent != -1) {
 447
 448                                 if (m->inUsersGroups(copy->tree[i].getName(), Groups)) {
 449                                         //find my siblings name
 450                                         int parentRC = copy->tree[parent].getRChild();
 451                                         int parentLC = copy->tree[parent].getLChild();
 452
 453                                         //if I am the right child, then my sib is the left child
 454                                         int sibIndex = parentRC;
 455                                         if (parentRC == i) { sibIndex = parentLC; }
 456
 457                                         string sibsName = copy->tree[sibIndex].getName();
 458
 459                                         //if yes, is my sibling
 460                                         if ((m->inUsersGroups(sibsName, Groups)) || (sibsName == "")) {
 461                                                 //we both are okay no trimming required
 462                                         }else{
 463                                                 //i am, my sib is not, so remove sib by setting my parent to my grandparent
 464                                                 int grandparent = copy->tree[parent].getParent();
 465                                                 int grandparentLC = copy->tree[grandparent].getLChild();
 466                                                 int grandparentRC = copy->tree[grandparent].getRChild();
 467
 468                                                 //whichever of my granparents children was my parent now equals me
 469                                                 if (grandparentLC == parent) { grandparentLC = i; }
 470                                                 else { grandparentRC = i; }
 471
 472                                                 copy->tree[i].setParent(grandparent);
 473                                                 copy->tree[i].setBranchLength((copy->tree[i].getBranchLength()+copy->tree[parent].getBranchLength()));
 474                                                 if (grandparent != -1) {
 475                                                         copy->tree[grandparent].setChildren(grandparentLC, grandparentRC);
 476                                                 }
 477                                                 removedLeaves.insert(sibIndex);
 478                                         }
 479                                 }else{
 480                                         //find my siblings name
 481                                         int parentRC = copy->tree[parent].getRChild();
 482                                         int parentLC = copy->tree[parent].getLChild();
 483
 484                                         //if I am the right child, then my sib is the left child
 485                                         int sibIndex = parentRC;
 486                                         if (parentRC == i) { sibIndex = parentLC; }
 487
 488                                         string sibsName = copy->tree[sibIndex].getName();
 489
 490                                         //if no is my sibling
 491                                         if ((m->inUsersGroups(sibsName, Groups)) || (sibsName == "")) {
 492                                                 //i am not, but my sib is
 493                                                 int grandparent = copy->tree[parent].getParent();
 494                                                 int grandparentLC = copy->tree[grandparent].getLChild();
 495                                                 int grandparentRC = copy->tree[grandparent].getRChild();
 496
 497                                                 //whichever of my granparents children was my parent now equals my sib
 498                                                 if (grandparentLC == parent) { grandparentLC = sibIndex; }
 499                                                 else { grandparentRC = sibIndex; }
 500
 501                                                 copy->tree[sibIndex].setParent(grandparent);
 502                                                 copy->tree[sibIndex].setBranchLength((copy->tree[sibIndex].getBranchLength()+copy->tree[parent].getBranchLength()));
 503                                                 if (grandparent != -1) {
 504                                                         copy->tree[grandparent].setChildren(grandparentLC, grandparentRC);
 505                                                 }
 506                                                 removedLeaves.insert(i);
 507                                         }else{
 508                                                 //neither of us are, so we want to eliminate ourselves and our parent
 509                                                 //so set our parents sib to our great-grandparent
 510                                                 int parent = copy->tree[i].getParent();
 511                                                 int grandparent = copy->tree[parent].getParent();
 512                                                 int parentsSibIndex;
 513                                                 if (grandparent != -1) {
 514                                                         int greatgrandparent = copy->tree[grandparent].getParent();
 515                                                         int greatgrandparentLC, greatgrandparentRC;
 516                                                         if (greatgrandparent != -1) {
 517                                                                 greatgrandparentLC = copy->tree[greatgrandparent].getLChild();
 518                                                                 greatgrandparentRC = copy->tree[greatgrandparent].getRChild();
 519                                                         }
 520
 521                                                         int grandparentLC = copy->tree[grandparent].getLChild();
 522                                                         int grandparentRC = copy->tree[grandparent].getRChild();
 523
 524                                                         parentsSibIndex = grandparentLC;
 525                                                         if (grandparentLC == parent) { parentsSibIndex = grandparentRC; }
 526
 527                                                         //whichever of my greatgrandparents children was my grandparent
 528                                                         if (greatgrandparentLC == grandparent) { greatgrandparentLC = parentsSibIndex; }
 529                                                         else { greatgrandparentRC = parentsSibIndex; }
 530
 531                                                         copy->tree[parentsSibIndex].setParent(greatgrandparent);
 532                                                         copy->tree[parentsSibIndex].setBranchLength((copy->tree[parentsSibIndex].getBranchLength()+copy->tree[grandparent].getBranchLength()));
 533                                                         if (greatgrandparent != -1) {
 534                                                                 copy->tree[greatgrandparent].setChildren(greatgrandparentLC, greatgrandparentRC);
 535                                                         }
 536                                                 }else{
 537                                                         copy->tree[parent].setParent(-1);
 538                                                         //cout << "issues with making subtree" << endl;
 539                                                 }
 540                                                 removedLeaves.insert(sibIndex);
 541                                                 removedLeaves.insert(i);
 542                                         }
 543                                 }
 544                         }
 545                         }
 546                 }
 547
 548                 int root = 0;
 549                 for (int i = 0; i < copy->getNumNodes(); i++) {
 550                         //you found the root
 551                         if (copy->tree[i].getParent() == -1) { root = i; break; }
 552                 }
 553
 554                 int nextSpot = numLeaves;
 555                 populateNewTree(copy->tree, root, nextSpot);
 556
 557         delete copy;
 558         }
 559         catch(exception& e) {
 560                 m->errorOut(e, "Tree", "getSubTree");
 561                 exit(1);
 562         }
 563 }
 564 /*****************************************************************
 565 //assumes nameMap contains unique names as key or is empty.
 566 //assumes numLeaves defined in tree constructor equals size of seqsToInclude and seqsToInclude only contains unique seqs.
 567 int Tree::getSubTree(Tree* copy, vector<string> seqsToInclude, map<string, string> nameMap) {
 568         try {
 569
 570         if (numLeaves != seqsToInclude.size()) { m->mothurOut("[ERROR]: numLeaves does not equal numUniques, cannot create subtree.\n"); m->control_pressed = true; return 0; }
 571
 572         getSubTree(copy, seqsToInclude);
 573         if (nameMap.size() != 0) {  addNamesToCounts(nameMap);  }
 574
 575         //build the pGroups in non leaf nodes to be used in the parsimony calcs.
 576                 for (int i = numLeaves; i < numNodes; i++) {
 577                         if (m->control_pressed) { return 1; }
 578
 579                         tree[i].pGroups = (mergeGroups(i));
 580                         tree[i].pcount = (mergeGcounts(i));
 581                 }
 582
 583         return 0;
 584     }
 585         catch(exception& e) {
 586                 m->errorOut(e, "Tree", "getSubTree");
 587                 exit(1);
 588         }
 589 }
 590 /*****************************************************************/
 591 int Tree::populateNewTree(vector<Node>& oldtree, int node, int& index) {
 592         try {
 593
 594                 if (oldtree[node].getLChild() != -1) {
 595                         int rc = populateNewTree(oldtree, oldtree[node].getLChild(), index);
 596                         int lc = populateNewTree(oldtree, oldtree[node].getRChild(), index);
 597
 598                         tree[index].setChildren(lc, rc);
 599                         tree[rc].setParent(index);
 600                         tree[lc].setParent(index);
 601
 602                         tree[index].setBranchLength(oldtree[node].getBranchLength());
 603                         tree[rc].setBranchLength(oldtree[oldtree[node].getLChild()].getBranchLength());
 604                         tree[lc].setBranchLength(oldtree[oldtree[node].getRChild()].getBranchLength());
 605
 606                         return (index++);
 607                 }else { //you are a leaf
 608                         int indexInNewTree = getIndex(oldtree[node].getName());
 609                         return indexInNewTree;
 610                 }
 611         }
 612         catch(exception& e) {
 613                 m->errorOut(e, "Tree", "populateNewTree");
 614                 exit(1);
 615         }
 616 }
 617 /*****************************************************************/
 618 void Tree::getCopy(Tree* copy, bool subsample) {
 619         try {
 620
 621                 //for each node in the tree copy its info
 622                 for (int i = 0; i < numNodes; i++) {
 623                         //copy branch length
 624                         tree[i].setBranchLength(copy->tree[i].getBranchLength());
 625
 626                         //copy parent
 627                         tree[i].setParent(copy->tree[i].getParent());
 628
 629                         //copy children
 630                         tree[i].setChildren(copy->tree[i].getLChild(), copy->tree[i].getRChild());
 631         }
 632
 633         //build the pGroups in non leaf nodes to be used in the parsimony calcs.
 634                 for (int i = numLeaves; i < numNodes; i++) {
 635                         if (m->control_pressed) { break; }
 636
 637                         tree[i].pGroups = (mergeGroups(i));
 638                         tree[i].pcount = (mergeGcounts(i));
 639                 }
 640         }
 641         catch(exception& e) {
 642                 m->errorOut(e, "Tree", "getCopy");
 643                 exit(1);
 644         }
 645 }
 646 /*****************************************************************/
 647 void Tree::getCopy(Tree* copy) {
 648         try {
 649
 650                 //for each node in the tree copy its info
 651                 for (int i = 0; i < numNodes; i++) {
 652                         //copy name
 653                         tree[i].setName(copy->tree[i].getName());
 654
 655                         //copy group
 656                         tree[i].setGroup(copy->tree[i].getGroup());
 657
 658                         //copy branch length
 659                         tree[i].setBranchLength(copy->tree[i].getBranchLength());
 660
 661                         //copy parent
 662                         tree[i].setParent(copy->tree[i].getParent());
 663
 664                         //copy children
 665                         tree[i].setChildren(copy->tree[i].getLChild(), copy->tree[i].getRChild());
 666
 667                         //copy index in node and tmap
 668             setIndex(copy->tree[i].getName(), getIndex(copy->tree[i].getName()));
 669                         tree[i].setIndex(copy->tree[i].getIndex());
 670
 671                         //copy pGroups
 672                         tree[i].pGroups = copy->tree[i].pGroups;
 673
 674                         //copy pcount
 675                         tree[i].pcount = copy->tree[i].pcount;
 676                 }
 677
 678                 groupNodeInfo = copy->groupNodeInfo;
 679
 680         }
 681         catch(exception& e) {
 682                 m->errorOut(e, "Tree", "getCopy");
 683                 exit(1);
 684         }
 685 }
 686 /*****************************************************************/
 687 //returns a map with a groupname and the number of times that group was seen in the children
 688 //for instance if your children are white and black then it would return a map with 2 entries
 689 // p[white] = 1 and p[black] = 1.  Now go up a level and merge that with a node who has p[white] = 1
 690 //and you get p[white] = 2, p[black] = 1, but you erase the p[black] because you have a p value higher than 1.
 691
 692 map<string, int> Tree::mergeGroups(int i) {
 693         try {
 694                 int lc = tree[i].getLChild();
 695                 int rc = tree[i].getRChild();
 696
 697                 //set parsimony groups to left child
 698                 map<string,int> parsimony = tree[lc].pGroups;
 699
 700                 int maxPars = 1;
 701
 702                 //look at right child groups and update maxPars if right child has something higher for that group.
 703                 for(it=tree[rc].pGroups.begin();it!=tree[rc].pGroups.end();it++){
 704                         it2 = parsimony.find(it->first);
 705                         if (it2 != parsimony.end()) {
 706                                 parsimony[it->first]++;
 707                         }else {
 708                                 parsimony[it->first] = 1;
 709                         }
 710
 711                         if(parsimony[it->first] > maxPars){
 712                                 maxPars = parsimony[it->first];
 713                         }
 714                 }
 715
 716                 // this is true if right child had a greater parsimony for a certain group
 717                 if(maxPars > 1){
 718                         //erase all the groups that are only 1 because you found something with 2.
 719                         for(it=parsimony.begin();it!=parsimony.end();){
 720                                 if(it->second == 1){
 721                                         parsimony.erase(it++);
 722                                 }else { it++; }
 723                         }
 724                         //set one remaining groups to 1
 725                         //so with our above example p[white] = 2 would be left and it would become p[white] = 1
 726                         for(it=parsimony.begin();it!=parsimony.end();it++){
 727                                 parsimony[it->first] = 1;
 728                         }
 729
 730                 }
 731
 732                 return parsimony;
 733         }
 734         catch(exception& e) {
 735                 m->errorOut(e, "Tree", "mergeGroups");
 736                 exit(1);
 737         }
 738 }
 739 /*****************************************************************/
 740 //returns a map with a groupname and the number of times that group was seen in the children
 741 //for instance if your children are white and black then it would return a map with 2 entries
 742 // p[white] = 1 and p[black] = 1.  Now go up a level and merge that with a node who has p[white] = 1
 743 //and you get p[white] = 2, p[black] = 1, but you erase the p[black] because you have a p value higher than 1.
 744
 745 map<string, int> Tree::mergeUserGroups(int i, vector<string> g) {
 746         try {
 747
 748                 int lc = tree[i].getLChild();
 749                 int rc = tree[i].getRChild();
 750
 751                 //loop through nodes groups removing the ones the user doesn't want
 752                 for(it=tree[lc].pGroups.begin();it!=tree[lc].pGroups.end();){
 753                                 if (m->inUsersGroups(it->first, g) != true) {
 754                                         tree[lc].pGroups.erase(it++);
 755                                 }else { it++; }
 756                 }
 757
 758                 //loop through nodes groups removing the ones the user doesn't want
 759                 for(it=tree[rc].pGroups.begin();it!=tree[rc].pGroups.end();){
 760                                 if (m->inUsersGroups(it->first, g) != true) {
 761                                         tree[rc].pGroups.erase(it++);
 762                                 }else { it++; }
 763                 }
 764
 765                 //set parsimony groups to left child
 766                 map<string,int> parsimony = tree[lc].pGroups;
 767
 768                 int maxPars = 1;
 769
 770                 //look at right child groups and update maxPars if right child has something higher for that group.
 771                 for(it=tree[rc].pGroups.begin();it!=tree[rc].pGroups.end();it++){
 772                         it2 = parsimony.find(it->first);
 773                         if (it2 != parsimony.end()) {
 774                                 parsimony[it->first]++;
 775                         }else {
 776                                 parsimony[it->first] = 1;
 777                         }
 778
 779                         if(parsimony[it->first] > maxPars){
 780                                 maxPars = parsimony[it->first];
 781                         }
 782                 }
 783
 784                 // this is true if right child had a greater parsimony for a certain group
 785                 if(maxPars > 1){
 786                         //erase all the groups that are only 1 because you found something with 2.
 787                         for(it=parsimony.begin();it!=parsimony.end();){
 788                                 if(it->second == 1){
 789                                         parsimony.erase(it++);
 790                                 }else { it++; }
 791                         }
 792
 793                         for(it=parsimony.begin();it!=parsimony.end();it++){
 794                                 parsimony[it->first] = 1;
 795                         }
 796                 }
 797
 798                 return parsimony;
 799         }
 800         catch(exception& e) {
 801                 m->errorOut(e, "Tree", "mergeUserGroups");
 802                 exit(1);
 803         }
 804 }
 805
 806
 807 /**************************************************************************************************/
 808
 809 map<string,int> Tree::mergeGcounts(int position) {
 810         try{
 811                 map<string,int>::iterator pos;
 812
 813                 int lc = tree[position].getLChild();
 814                 int rc = tree[position].getRChild();
 815
 816                 map<string,int> sum = tree[lc].pcount;
 817
 818                 for(it=tree[rc].pcount.begin();it!=tree[rc].pcount.end();it++){
 819                         sum[it->first] += it->second;
 820                 }
 821                 return sum;
 822         }
 823         catch(exception& e) {
 824                 m->errorOut(e, "Tree", "mergeGcounts");
 825                 exit(1);
 826         }
 827 }
 828 /**************************************************************************************************/
 829 void Tree::randomLabels(vector<string> g) {
 830         try {
 831
 832                 //initialize groupNodeInfo
 833                 for (int i = 0; i < (ct->getNamesOfGroups()).size(); i++) {
 834                         groupNodeInfo[(ct->getNamesOfGroups())[i]].resize(0);
 835                 }
 836
 837                 for(int i = 0; i < numLeaves; i++){
 838                         int z;
 839                         //get random index to switch with
 840                         z = int((float)(i+1) * (float)(rand()) / ((float)RAND_MAX+1.0));
 841
 842                         //you only want to randomize the nodes that are from a group the user wants analyzed, so
 843                         //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.
 844                         bool treez, treei;
 845
 846                         treez = m->inUsersGroups(tree[z].getGroup(), g);
 847                         treei = m->inUsersGroups(tree[i].getGroup(), g);
 848
 849                         if ((treez == true) && (treei == true)) {
 850                                 //switches node i and node z's info.
 851                                 map<string,int> lib_hold = tree[z].pGroups;
 852                                 tree[z].pGroups = (tree[i].pGroups);
 853                                 tree[i].pGroups = (lib_hold);
 854
 855                                 vector<string> zgroup = tree[z].getGroup();
 856                                 tree[z].setGroup(tree[i].getGroup());
 857                                 tree[i].setGroup(zgroup);
 858
 859                                 string zname = tree[z].getName();
 860                                 tree[z].setName(tree[i].getName());
 861                                 tree[i].setName(zname);
 862
 863                                 map<string,int> gcount_hold = tree[z].pcount;
 864                                 tree[z].pcount = (tree[i].pcount);
 865                                 tree[i].pcount = (gcount_hold);
 866                         }
 867
 868                         for (int k = 0; k < (tree[i].getGroup()).size(); k++) {  groupNodeInfo[(tree[i].getGroup())[k]].push_back(i); }
 869                         for (int k = 0; k < (tree[z].getGroup()).size(); k++) {  groupNodeInfo[(tree[z].getGroup())[k]].push_back(z); }
 870                 }
 871         }
 872         catch(exception& e) {
 873                 m->errorOut(e, "Tree", "randomLabels");
 874                 exit(1);
 875         }
 876 }
 877 /**************************************************************************************************/
 878 void Tree::randomBlengths()  {
 879         try {
 880                 for(int i=numNodes-1;i>=0;i--){
 881                         int z = int((float)(i+1) * (float)(rand()) / ((float)RAND_MAX+1.0));
 882
 883                         float bl_hold = tree[z].getBranchLength();
 884                         tree[z].setBranchLength(tree[i].getBranchLength());
 885                         tree[i].setBranchLength(bl_hold);
 886                 }
 887         }
 888         catch(exception& e) {
 889                 m->errorOut(e, "Tree", "randomBlengths");
 890                 exit(1);
 891         }
 892 }
 893 /*************************************************************************************************/
 894 void Tree::assembleRandomUnifracTree(vector<string> g) {
 895         randomLabels(g);
 896         assembleTree();
 897 }
 898 /*************************************************************************************************/
 899 void Tree::assembleRandomUnifracTree(string groupA, string groupB) {
 900         vector<string> temp; temp.push_back(groupA); temp.push_back(groupB);
 901         randomLabels(temp);
 902         assembleTree();
 903 }
 904
 905 /*************************************************************************************************/
 906 //for now it's just random topology but may become random labels as well later that why this is such a simple function now...
 907 void Tree::assembleRandomTree() {
 908         randomTopology();
 909         assembleTree();
 910 }
 911 /**************************************************************************************************/
 912
 913 void Tree::randomTopology() {
 914         try {
 915                 for(int i=0;i<numNodes;i++){
 916                         tree[i].setParent(-1);
 917                 }
 918                 for(int i=numLeaves;i<numNodes;i++){
 919                         tree[i].setChildren(-1, -1);
 920                 }
 921
 922                 for(int i=numLeaves;i<numNodes;i++){
 923                         int escape =0;
 924                         int rnd_index1, rnd_index2;
 925                         while(escape == 0){
 926                                 rnd_index1 = (int)(((double)rand() / (double) RAND_MAX)*i);
 927                                 if(tree[rnd_index1].getParent() == -1){escape = 1;}
 928                         }
 929
 930                         escape = 0;
 931                         while(escape == 0){
 932                                 rnd_index2 = (int)(((double)rand() / (double) RAND_MAX)*i);
 933                                 if(rnd_index2 != rnd_index1 && tree[rnd_index2].getParent() == -1){
 934                                         escape = 1;
 935                                 }
 936                         }
 937
 938                         tree[i].setChildren(rnd_index1,rnd_index2);
 939                         tree[i].setParent(-1);
 940                         tree[rnd_index1].setParent(i);
 941                         tree[rnd_index2].setParent(i);
 942                 }
 943         }
 944         catch(exception& e) {
 945                 m->errorOut(e, "Tree", "randomTopology");
 946                 exit(1);
 947         }
 948 }
 949 /*****************************************************************/
 950 void Tree::print(ostream& out) {
 951         try {
 952                 int root = findRoot();
 953                 printBranch(root, out, "branch");
 954                 out << ";" << endl;
 955         }
 956         catch(exception& e) {
 957                 m->errorOut(e, "Tree", "print");
 958                 exit(1);
 959         }
 960 }
 961 /*****************************************************************/
 962 void Tree::print(ostream& out, map<string, string> nameMap) {
 963         try {
 964                 int root = findRoot();
 965                 printBranch(root, out, nameMap);
 966                 out << ";" << endl;
 967         }
 968         catch(exception& e) {
 969                 m->errorOut(e, "Tree", "print");
 970                 exit(1);
 971         }
 972 }
 973 /*****************************************************************/
 974 void Tree::print(ostream& out, string mode) {
 975         try {
 976                 int root = findRoot();
 977                 printBranch(root, out, mode);
 978                 out << ";" << endl;
 979         }
 980         catch(exception& e) {
 981                 m->errorOut(e, "Tree", "print");
 982                 exit(1);
 983         }
 984 }
 985 /*****************************************************************/
 986 // This prints out the tree in Newick form.
 987 void Tree::createNewickFile(string f) {
 988         try {
 989                 int root = findRoot();
 990
 991                 filename = f;
 992
 993                 m->openOutputFile(filename, out);
 994
 995                 printBranch(root, out, "branch");
 996
 997                 // you are at the end of the tree
 998                 out << ";" << endl;
 999                 out.close();
1000         }
1001         catch(exception& e) {
1002                 m->errorOut(e, "Tree", "createNewickFile");
1003                 exit(1);
1004         }
1005 }
1006
1007 /*****************************************************************/
1008 //This function finds the index of the root node.
1009
1010 int Tree::findRoot() {
1011         try {
1012                 for (int i = 0; i < numNodes; i++) {
1013                         //you found the root
1014                         if (tree[i].getParent() == -1) { return i; }
1015                         //cout << "i = " << i << endl;
1016                         //cout << "i's parent = " << tree[i].getParent() << endl;
1017                 }
1018                 return -1;
1019         }
1020         catch(exception& e) {
1021                 m->errorOut(e, "Tree", "findRoot");
1022                 exit(1);
1023         }
1024 }
1025 /*****************************************************************/
1026 void Tree::printBranch(int node, ostream& out, map<string, string> names) {
1027 try {
1028
1029 // you are not a leaf
1030                 if (tree[node].getLChild() != -1) {
1031                         out << "(";
1032                         printBranch(tree[node].getLChild(), out, names);
1033                         out << ",";
1034                         printBranch(tree[node].getRChild(), out, names);
1035                         out << ")";
1036
1037             //if there is a branch length then print it
1038             if (tree[node].getBranchLength() != -1) {
1039                 out << ":" << tree[node].getBranchLength();
1040             }
1041
1042                 }else { //you are a leaf
1043             map<string, string>::iterator itNames = names.find(tree[node].getName());
1044
1045             string outputString = "";
1046             if (itNames != names.end()) {
1047
1048                 vector<string> dupNames;
1049                 m->splitAtComma((itNames->second), dupNames);
1050
1051                 if (dupNames.size() == 1) {
1052                     outputString += tree[node].getName();
1053                     if (tree[node].getBranchLength() != -1) {
1054                         outputString += ":" + toString(tree[node].getBranchLength());
1055                     }
1056                 }else {
1057                     outputString += "(";
1058
1059                     for (int u = 0; u < dupNames.size()-1; u++) {
1060                         outputString += dupNames[u];
1061
1062                         if (tree[node].getBranchLength() != -1) {
1063                             outputString += ":" + toString(0.0);
1064                         }
1065                         outputString += ",";
1066                     }
1067
1068                     outputString += dupNames[dupNames.size()-1];
1069                     if (tree[node].getBranchLength() != -1) {
1070                         outputString += ":" + toString(0.0);
1071                     }
1072
1073                     outputString += ")";
1074                     if (tree[node].getBranchLength() != -1) {
1075                         outputString += ":" + toString(tree[node].getBranchLength());
1076                     }
1077                 }
1078             }else {
1079                 outputString = tree[node].getName();
1080                 //if there is a branch length then print it
1081                 if (tree[node].getBranchLength() != -1) {
1082                     outputString += ":" + toString(tree[node].getBranchLength());
1083                 }
1084
1085                 m->mothurOut("[ERROR]: " + tree[node].getName() + " is not in your namefile, please correct."); m->mothurOutEndLine();
1086             }
1087
1088             out << outputString;
1089                 }
1090
1091         }
1092         catch(exception& e) {
1093                 m->errorOut(e, "Tree", "printBranch");
1094                 exit(1);
1095         }
1096 }
1097 /*****************************************************************/
1098 void Tree::printBranch(int node, ostream& out, string mode) {
1099     try {
1100
1101         // you are not a leaf
1102                 if (tree[node].getLChild() != -1) {
1103                         out << "(";
1104                         printBranch(tree[node].getLChild(), out, mode);
1105                         out << ",";
1106                         printBranch(tree[node].getRChild(), out, mode);
1107                         out << ")";
1108                         if (mode == "branch") {
1109                                 //if there is a branch length then print it
1110                                 if (tree[node].getBranchLength() != -1) {
1111                                         out << ":" << tree[node].getBranchLength();
1112                                 }
1113                         }else if (mode == "boot") {
1114                                 //if there is a label then print it
1115                                 if (tree[node].getLabel() != -1) {
1116                                         out << tree[node].getLabel();
1117                                 }
1118                         }else if (mode == "both") {
1119                                 if (tree[node].getLabel() != -1) {
1120                                         out << tree[node].getLabel();
1121                                 }
1122                                 //if there is a branch length then print it
1123                                 if (tree[node].getBranchLength() != -1) {
1124                                         out << ":" << tree[node].getBranchLength();
1125                                 }
1126                         }
1127                 }else { //you are a leaf
1128                         vector<string> leafGroup = ct->getGroups(tree[node].getName());
1129
1130                         if (mode == "branch") {
1131                                 out << leafGroup[0];
1132                                 //if there is a branch length then print it
1133                                 if (tree[node].getBranchLength() != -1) {
1134                                         out << ":" << tree[node].getBranchLength();
1135                                 }
1136                         }else if (mode == "boot") {
1137                                 out << leafGroup[0];
1138                                 //if there is a label then print it
1139                                 if (tree[node].getLabel() != -1) {
1140                                         out << tree[node].getLabel();
1141                                 }
1142                         }else if (mode == "both") {
1143                                 out << tree[node].getName();
1144                                 if (tree[node].getLabel() != -1) {
1145                                         out << tree[node].getLabel();
1146                                 }
1147                                 //if there is a branch length then print it
1148                                 if (tree[node].getBranchLength() != -1) {
1149                                         out << ":" << tree[node].getBranchLength();
1150                                 }
1151                         }
1152                 }
1153
1154         }
1155         catch(exception& e) {
1156                 m->errorOut(e, "Tree", "printBranch");
1157                 exit(1);
1158         }
1159 }
1160 /*****************************************************************/
1161 void Tree::printBranch(int node, ostream& out, string mode, vector<Node>& theseNodes) {
1162         try {
1163
1164                 // you are not a leaf
1165                 if (theseNodes[node].getLChild() != -1) {
1166                         out << "(";
1167                         printBranch(theseNodes[node].getLChild(), out, mode);
1168                         out << ",";
1169                         printBranch(theseNodes[node].getRChild(), out, mode);
1170                         out << ")";
1171                         if (mode == "branch") {
1172                                 //if there is a branch length then print it
1173                                 if (theseNodes[node].getBranchLength() != -1) {
1174                                         out << ":" << theseNodes[node].getBranchLength();
1175                                 }
1176                         }else if (mode == "boot") {
1177                                 //if there is a label then print it
1178                                 if (theseNodes[node].getLabel() != -1) {
1179                                         out << theseNodes[node].getLabel();
1180                                 }
1181                         }else if (mode == "both") {
1182                                 if (theseNodes[node].getLabel() != -1) {
1183                                         out << theseNodes[node].getLabel();
1184                                 }
1185                                 //if there is a branch length then print it
1186                                 if (theseNodes[node].getBranchLength() != -1) {
1187                                         out << ":" << theseNodes[node].getBranchLength();
1188                                 }
1189                         }
1190                 }else { //you are a leaf
1191                         vector<string> leafGroup = ct->getGroups(theseNodes[node].getName());
1192
1193                         if (mode == "branch") {
1194                                 out << leafGroup[0];
1195                                 //if there is a branch length then print it
1196                                 if (theseNodes[node].getBranchLength() != -1) {
1197                                         out << ":" << theseNodes[node].getBranchLength();
1198                                 }
1199                         }else if (mode == "boot") {
1200                                 out << leafGroup[0];
1201                                 //if there is a label then print it
1202                                 if (theseNodes[node].getLabel() != -1) {
1203                                         out << theseNodes[node].getLabel();
1204                                 }
1205                         }else if (mode == "both") {
1206                                 out << theseNodes[node].getName();
1207                                 if (theseNodes[node].getLabel() != -1) {
1208                                         out << theseNodes[node].getLabel();
1209                                 }
1210                                 //if there is a branch length then print it
1211                                 if (theseNodes[node].getBranchLength() != -1) {
1212                                         out << ":" << theseNodes[node].getBranchLength();
1213                                 }
1214                         }
1215                 }
1216
1217         }
1218         catch(exception& e) {
1219                 m->errorOut(e, "Tree", "printBranch");
1220                 exit(1);
1221         }
1222 }
1223 /*****************************************************************/
1224
1225 void Tree::printTree() {
1226
1227         for(int i=0;i<numNodes;i++){
1228                 cout << i << '\t';
1229                 tree[i].printNode();
1230         }
1231
1232 }
1233
1234 /*****************************************************************/
1235 //this code is a mess and should be rethought...-slw
1236 int Tree::parseTreeFile() {
1237
1238         //only takes names from the first tree and assumes that all trees use the same names.
1239         try {
1240                 string filename = m->getTreeFile();
1241                 ifstream filehandle;
1242                 m->openInputFile(filename, filehandle);
1243                 int c, comment;
1244                 comment = 0;
1245                 int done = 1;
1246
1247                 //ifyou are not a nexus file
1248                 if((c = filehandle.peek()) != '#') {
1249                         while((c = filehandle.peek()) != ';') {
1250                 if (m->control_pressed) {  filehandle.close(); return 0; }
1251                                 while ((c = filehandle.peek()) != ';') {
1252                     if (m->control_pressed) {  filehandle.close(); return 0; }
1253                                         // get past comments
1254                                         if(c == '[') {
1255                                                 comment = 1;
1256                                         }
1257                                         if(c == ']'){
1258                                                 comment = 0;
1259                                         }
1260                                         if((c == '(') && (comment != 1)){ break; }
1261                                         filehandle.get();
1262                                 }
1263
1264                                 done = readTreeString(filehandle);
1265                                 if (done == 0) { break; }
1266                         }
1267                 //ifyou are a nexus file
1268                 }else if((c = filehandle.peek()) == '#') {
1269                         string holder = "";
1270
1271                         // get past comments
1272                         while(holder != "translate" && holder != "Translate"){
1273                 if (m->control_pressed) {  filehandle.close(); return 0; }
1274                                 if(holder == "[" || holder == "[!"){
1275                                         comment = 1;
1276                                 }
1277                                 if(holder == "]"){
1278                                         comment = 0;
1279                                 }
1280                                 filehandle >> holder;
1281
1282                                 //if there is no translate then you must read tree string otherwise use translate to get names
1283                                 if((holder == "tree") && (comment != 1)){
1284                                         //pass over the "tree rep.6878900 = "
1285                                         while (((c = filehandle.get()) != '(') && ((c = filehandle.peek()) != EOF)) {;}
1286
1287                                         if(c == EOF) { break; }
1288                                         filehandle.putback(c);  //put back first ( of tree.
1289                                         done = readTreeString(filehandle);
1290
1291                                         break;
1292                                 }
1293
1294                                 if (done == 0) { break;  }
1295                         }
1296
1297                         //use nexus translation rather than parsing tree to save time
1298                         if((holder == "translate") || (holder == "Translate")) {
1299
1300                                 string number, name, h;
1301                                 h = ""; // so it enters the loop the first time
1302                                 while((h != ";") && (number != ";")) {
1303                     if (m->control_pressed) {  filehandle.close(); return 0; }
1304                                         filehandle >> number;
1305                                         filehandle >> name;
1306
1307                                         //c = , until done with translation then c = ;
1308                                         h = name.substr(name.length()-1, name.length());
1309                                         name.erase(name.end()-1);  //erase the comma
1310                                         m->Treenames.push_back(number);
1311                                 }
1312                                 if(number == ";") { m->Treenames.pop_back(); }  //in case ';' from translation is on next line instead of next to last name
1313                         }
1314                 }
1315                 filehandle.close();
1316                 return 0;
1317                 //for (int i = 0; i < globaldata->Treenames.size(); i++) {
1318 //cout << globaldata->Treenames[i] << endl; }
1319 //cout << globaldata->Treenames.size() << endl;
1320         }
1321         catch(exception& e) {
1322                 m->errorOut(e, "Tree", "parseTreeFile");
1323                 exit(1);
1324         }
1325 }
1326 /*******************************************************/
1327
1328 /*******************************************************/
1329 int Tree::readTreeString(ifstream& filehandle)  {
1330         try {
1331                 int c;
1332                 string name;  //, k
1333
1334                 while((c = filehandle.peek()) != ';') {
1335             if (m->control_pressed) {  return 0; }
1336 //k = c;
1337 //cout << " at beginning of while " <<  k << endl;
1338                         if(c == ')')  {
1339                                 //to pass over labels in trees
1340                                 c=filehandle.get();
1341                                 while((c!=',') && (c != -1) && (c!= ':') && (c!=';')){ c=filehandle.get(); }
1342                                 filehandle.putback(c);
1343                         }
1344                         if(c == ';') { return 0; }
1345                         if(c == -1) { return 0; }
1346                         //if you are a name
1347                         if((c != '(') && (c != ')') && (c != ',') && (c != ':') && (c != '\n') && (c != '\t') && (c != 32)) { //32 is space
1348                                 name = "";
1349                                 c = filehandle.get();
1350                         //k = c;
1351 //cout << k << endl;
1352                                 while ((c != '(') && (c != ')') && (c != ',') && (c != ':')  && (c != '\n') && (c != 32) && (c != '\t')) {
1353                                         name += c;
1354                                         c = filehandle.get();
1355                         //k = c;
1356 //cout << " in name while " << k << endl;
1357                                 }
1358
1359 //cout << "name = " << name << endl;
1360                 if (name != "\r" ) {
1361                     m->Treenames.push_back(name);   } //cout << m->Treenames.size() << '\t' << name << endl;
1362
1363                                 filehandle.putback(c);
1364 //k = c;
1365 //cout << " after putback" <<  k << endl;
1366                         }
1367
1368                         if(c  == ':') { //read until you reach the end of the branch length
1369                                 while ((c != '(') && (c != ')') && (c != ',') && (c != ';') && (c != '\n') && (c != '\t') && (c != 32)) {
1370                                         c = filehandle.get();
1371         //k = c;
1372         //cout << " in branch while " << k << endl;
1373                                 }
1374                                 filehandle.putback(c);
1375                         }
1376
1377                         c = filehandle.get();
1378 //k = c;
1379         //cout << " here after get " << k << endl;
1380                         if(c == ';') { return 0; }
1381                         if(c == ')') { filehandle.putback(c); }
1382         //k = c;
1383 //cout << k << endl;
1384
1385                 }
1386                 return 0;
1387         }
1388         catch(exception& e) {
1389                 m->errorOut(e, "Tree", "readTreeString");
1390                 exit(1);
1391         }
1392 }
1393
1394 /*******************************************************/
1395
1396 /*******************************************************/
1397