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