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