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