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