]> git.donarmstrong.com Git - mothur.git/blob - tree.cpp
Revert to previous commit
[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) {
592         try {
593         
594                 //for each node in the tree copy its info
595                 for (int i = 0; i < numNodes; i++) {
596                         //copy name
597                         tree[i].setName(copy->tree[i].getName());
598                 
599                         //copy group
600                         tree[i].setGroup(copy->tree[i].getGroup());
601                         
602                         //copy branch length
603                         tree[i].setBranchLength(copy->tree[i].getBranchLength());
604                 
605                         //copy parent
606                         tree[i].setParent(copy->tree[i].getParent());
607                 
608                         //copy children
609                         tree[i].setChildren(copy->tree[i].getLChild(), copy->tree[i].getRChild());
610                 
611                         //copy index in node and tmap
612                         tree[i].setIndex(copy->tree[i].getIndex());
613                         setIndex(copy->tree[i].getName(), getIndex(copy->tree[i].getName()));
614                         
615                         //copy pGroups
616                         tree[i].pGroups = copy->tree[i].pGroups;
617                 
618                         //copy pcount
619                         tree[i].pcount = copy->tree[i].pcount;
620                 }
621                 
622                 groupNodeInfo = copy->groupNodeInfo;
623                 
624         }
625         catch(exception& e) {
626                 m->errorOut(e, "Tree", "getCopy");
627                 exit(1);
628         }
629 }
630 /*****************************************************************/
631 //returns a map with a groupname and the number of times that group was seen in the children
632 //for instance if your children are white and black then it would return a map with 2 entries
633 // p[white] = 1 and p[black] = 1.  Now go up a level and merge that with a node who has p[white] = 1
634 //and you get p[white] = 2, p[black] = 1, but you erase the p[black] because you have a p value higher than 1.
635
636 map<string, int> Tree::mergeGroups(int i) {
637         try {
638                 int lc = tree[i].getLChild();
639                 int rc = tree[i].getRChild();
640
641                 //set parsimony groups to left child
642                 map<string,int> parsimony = tree[lc].pGroups;
643                 
644                 int maxPars = 1;
645
646                 //look at right child groups and update maxPars if right child has something higher for that group.
647                 for(it=tree[rc].pGroups.begin();it!=tree[rc].pGroups.end();it++){
648                         it2 = parsimony.find(it->first);
649                         if (it2 != parsimony.end()) {
650                                 parsimony[it->first]++;
651                         }else {
652                                 parsimony[it->first] = 1;
653                         }
654                         
655                         if(parsimony[it->first] > maxPars){
656                                 maxPars = parsimony[it->first];
657                         }
658                 }
659         
660                 // this is true if right child had a greater parsimony for a certain group
661                 if(maxPars > 1){
662                         //erase all the groups that are only 1 because you found something with 2.
663                         for(it=parsimony.begin();it!=parsimony.end();){
664                                 if(it->second == 1){
665                                         parsimony.erase(it++);
666                                 }else { it++; }
667                         }
668                         //set one remaining groups to 1
669                         //so with our above example p[white] = 2 would be left and it would become p[white] = 1
670                         for(it=parsimony.begin();it!=parsimony.end();it++){
671                                 parsimony[it->first] = 1;
672                         }
673                 
674                 }
675         
676                 return parsimony;
677         }
678         catch(exception& e) {
679                 m->errorOut(e, "Tree", "mergeGroups");
680                 exit(1);
681         }
682 }
683 /*****************************************************************/
684 //returns a map with a groupname and the number of times that group was seen in the children
685 //for instance if your children are white and black then it would return a map with 2 entries
686 // p[white] = 1 and p[black] = 1.  Now go up a level and merge that with a node who has p[white] = 1
687 //and you get p[white] = 2, p[black] = 1, but you erase the p[black] because you have a p value higher than 1.
688
689 map<string, int> Tree::mergeUserGroups(int i, vector<string> g) {
690         try {
691         
692                 int lc = tree[i].getLChild();
693                 int rc = tree[i].getRChild();
694                 
695                 //loop through nodes groups removing the ones the user doesn't want
696                 for(it=tree[lc].pGroups.begin();it!=tree[lc].pGroups.end();){
697                                 if (m->inUsersGroups(it->first, g) != true) {
698                                         tree[lc].pGroups.erase(it++);
699                                 }else { it++; }
700                 }
701
702                 //loop through nodes groups removing the ones the user doesn't want
703                 for(it=tree[rc].pGroups.begin();it!=tree[rc].pGroups.end();){
704                                 if (m->inUsersGroups(it->first, g) != true) {
705                                         tree[rc].pGroups.erase(it++);
706                                 }else { it++; }
707                 }
708
709                 //set parsimony groups to left child
710                 map<string,int> parsimony = tree[lc].pGroups;
711                 
712                 int maxPars = 1;
713
714                 //look at right child groups and update maxPars if right child has something higher for that group.
715                 for(it=tree[rc].pGroups.begin();it!=tree[rc].pGroups.end();it++){
716                         it2 = parsimony.find(it->first);
717                         if (it2 != parsimony.end()) {
718                                 parsimony[it->first]++;
719                         }else {
720                                 parsimony[it->first] = 1;
721                         }
722                         
723                         if(parsimony[it->first] > maxPars){
724                                 maxPars = parsimony[it->first];
725                         }
726                 }
727                         
728                 // this is true if right child had a greater parsimony for a certain group
729                 if(maxPars > 1){
730                         //erase all the groups that are only 1 because you found something with 2.
731                         for(it=parsimony.begin();it!=parsimony.end();){
732                                 if(it->second == 1){
733                                         parsimony.erase(it++);
734                                 }else { it++; }
735                         }
736
737                         for(it=parsimony.begin();it!=parsimony.end();it++){
738                                 parsimony[it->first] = 1;
739                         }
740                 }               
741                 
742                 return parsimony;
743         }
744         catch(exception& e) {
745                 m->errorOut(e, "Tree", "mergeUserGroups");
746                 exit(1);
747         }
748 }
749
750
751 /**************************************************************************************************/
752
753 map<string,int> Tree::mergeGcounts(int position) {
754         try{
755                 map<string,int>::iterator pos;
756         
757                 int lc = tree[position].getLChild();
758                 int rc = tree[position].getRChild();
759         
760                 map<string,int> sum = tree[lc].pcount;
761     
762                 for(it=tree[rc].pcount.begin();it!=tree[rc].pcount.end();it++){
763                         sum[it->first] += it->second;
764                 }
765                 return sum;
766         }
767         catch(exception& e) {
768                 m->errorOut(e, "Tree", "mergeGcounts");
769                 exit(1);
770         }
771 }
772 /**************************************************************************************************/
773 void Tree::randomLabels(vector<string> g) {
774         try {
775         
776                 //initialize groupNodeInfo
777                 for (int i = 0; i < (tmap->getNamesOfGroups()).size(); i++) {
778                         groupNodeInfo[(tmap->getNamesOfGroups())[i]].resize(0);
779                 }
780                 
781                 for(int i = 0; i < numLeaves; i++){
782                         int z;
783                         //get random index to switch with
784                         z = int((float)(i+1) * (float)(rand()) / ((float)RAND_MAX+1.0));        
785                         
786                         //you only want to randomize the nodes that are from a group the user wants analyzed, so
787                         //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.
788                         bool treez, treei;
789                 
790                         treez = m->inUsersGroups(tree[z].getGroup(), g);
791                         treei = m->inUsersGroups(tree[i].getGroup(), g);
792                         
793                         if ((treez == true) && (treei == true)) {
794                                 //switches node i and node z's info.
795                                 map<string,int> lib_hold = tree[z].pGroups;
796                                 tree[z].pGroups = (tree[i].pGroups);
797                                 tree[i].pGroups = (lib_hold);
798                                 
799                                 vector<string> zgroup = tree[z].getGroup();
800                                 tree[z].setGroup(tree[i].getGroup());
801                                 tree[i].setGroup(zgroup);
802                                 
803                                 string zname = tree[z].getName();
804                                 tree[z].setName(tree[i].getName());
805                                 tree[i].setName(zname);
806                                 
807                                 map<string,int> gcount_hold = tree[z].pcount;
808                                 tree[z].pcount = (tree[i].pcount);
809                                 tree[i].pcount = (gcount_hold);
810                         }
811                         
812                         for (int k = 0; k < (tree[i].getGroup()).size(); k++) {  groupNodeInfo[(tree[i].getGroup())[k]].push_back(i); }
813                         for (int k = 0; k < (tree[z].getGroup()).size(); k++) {  groupNodeInfo[(tree[z].getGroup())[k]].push_back(z); }
814                 }
815         }
816         catch(exception& e) {
817                 m->errorOut(e, "Tree", "randomLabels");
818                 exit(1);
819         }
820 }
821 /**************************************************************************************************/
822 void Tree::randomBlengths()  {
823         try {
824                 for(int i=numNodes-1;i>=0;i--){
825                         int z = int((float)(i+1) * (float)(rand()) / ((float)RAND_MAX+1.0));    
826                 
827                         float bl_hold = tree[z].getBranchLength();
828                         tree[z].setBranchLength(tree[i].getBranchLength());
829                         tree[i].setBranchLength(bl_hold);
830                 }
831         }
832         catch(exception& e) {
833                 m->errorOut(e, "Tree", "randomBlengths");
834                 exit(1);
835         }
836 }
837 /*************************************************************************************************/
838 void Tree::assembleRandomUnifracTree(vector<string> g) {
839         randomLabels(g);
840     map<string, string> empty;
841         assembleTree(empty);
842 }
843 /*************************************************************************************************/
844 void Tree::assembleRandomUnifracTree(string groupA, string groupB) {
845         vector<string> temp; temp.push_back(groupA); temp.push_back(groupB);
846         randomLabels(temp);
847     map<string, string> empty;
848         assembleTree(empty);
849 }
850
851 /*************************************************************************************************/
852 //for now it's just random topology but may become random labels as well later that why this is such a simple function now...
853 void Tree::assembleRandomTree() {
854         randomTopology();
855     map<string, string> empty;
856         assembleTree(empty);
857 }
858 /**************************************************************************************************/
859
860 void Tree::randomTopology() {
861         try {
862                 for(int i=0;i<numNodes;i++){
863                         tree[i].setParent(-1);
864                 }
865                 for(int i=numLeaves;i<numNodes;i++){
866                         tree[i].setChildren(-1, -1); 
867                 }
868     
869                 for(int i=numLeaves;i<numNodes;i++){
870                         int escape =0;
871                         int rnd_index1, rnd_index2;
872                         while(escape == 0){
873                                 rnd_index1 = (int)(((double)rand() / (double) RAND_MAX)*i);
874                                 if(tree[rnd_index1].getParent() == -1){escape = 1;}
875                         }
876                 
877                         escape = 0;
878                         while(escape == 0){
879                                 rnd_index2 = (int)(((double)rand() / (double) RAND_MAX)*i);
880                                 if(rnd_index2 != rnd_index1 && tree[rnd_index2].getParent() == -1){
881                                         escape = 1;
882                                 }               
883                         }
884         
885                         tree[i].setChildren(rnd_index1,rnd_index2);
886                         tree[i].setParent(-1);
887                         tree[rnd_index1].setParent(i);
888                         tree[rnd_index2].setParent(i);
889                 }
890         }
891         catch(exception& e) {
892                 m->errorOut(e, "Tree", "randomTopology");
893                 exit(1);
894         }
895 }
896 /*****************************************************************/
897 void Tree::print(ostream& out) {
898         try {
899                 int root = findRoot();
900                 printBranch(root, out, "branch");
901                 out << ";" << endl;
902         }
903         catch(exception& e) {
904                 m->errorOut(e, "Tree", "print");
905                 exit(1);
906         }
907 }
908 /*****************************************************************/
909 void Tree::print(ostream& out, map<string, string> nameMap) {
910         try {
911                 int root = findRoot();
912                 printBranch(root, out, nameMap);
913                 out << ";" << endl;
914         }
915         catch(exception& e) {
916                 m->errorOut(e, "Tree", "print");
917                 exit(1);
918         }
919 }
920 /*****************************************************************/
921 void Tree::print(ostream& out, string mode) {
922         try {
923                 int root = findRoot();
924                 printBranch(root, out, mode);
925                 out << ";" << endl;
926         }
927         catch(exception& e) {
928                 m->errorOut(e, "Tree", "print");
929                 exit(1);
930         }
931 }
932 /*****************************************************************/
933 // This prints out the tree in Newick form.
934 void Tree::createNewickFile(string f) {
935         try {
936                 int root = findRoot();
937         
938                 filename = f;
939
940                 m->openOutputFile(filename, out);
941                 
942                 printBranch(root, out, "branch");
943                 
944                 // you are at the end of the tree
945                 out << ";" << endl;
946                 out.close();
947         }
948         catch(exception& e) {
949                 m->errorOut(e, "Tree", "createNewickFile");
950                 exit(1);
951         }
952 }
953
954 /*****************************************************************/
955 //This function finds the index of the root node.
956
957 int Tree::findRoot() {
958         try {
959                 for (int i = 0; i < numNodes; i++) {
960                         //you found the root
961                         if (tree[i].getParent() == -1) { return i; }
962                         //cout << "i = " << i << endl;
963                         //cout << "i's parent = " << tree[i].getParent() << endl;  
964                 }
965                 return -1;
966         }
967         catch(exception& e) {
968                 m->errorOut(e, "Tree", "findRoot");
969                 exit(1);
970         }
971 }
972 /*****************************************************************/
973 void Tree::printBranch(int node, ostream& out, map<string, string> names) {
974 try {
975
976 // you are not a leaf
977                 if (tree[node].getLChild() != -1) {
978                         out << "(";
979                         printBranch(tree[node].getLChild(), out, names);
980                         out << ",";
981                         printBranch(tree[node].getRChild(), out, names);
982                         out << ")";
983                         
984             //if there is a branch length then print it
985             if (tree[node].getBranchLength() != -1) {
986                 out << ":" << tree[node].getBranchLength();
987             }
988                         
989                 }else { //you are a leaf
990             map<string, string>::iterator itNames = names.find(tree[node].getName());
991             
992             string outputString = "";
993             if (itNames != names.end()) { 
994                 
995                 vector<string> dupNames;
996                 m->splitAtComma((itNames->second), dupNames);
997                 
998                 if (dupNames.size() == 1) {
999                     outputString += tree[node].getName();
1000                     if (tree[node].getBranchLength() != -1) {
1001                         outputString += ":" + toString(tree[node].getBranchLength());
1002                     }
1003                 }else {
1004                     outputString += "(";
1005                     
1006                     for (int u = 0; u < dupNames.size()-1; u++) {
1007                         outputString += dupNames[u];
1008                         
1009                         if (tree[node].getBranchLength() != -1) {
1010                             outputString += ":" + toString(0.0);
1011                         }
1012                         outputString += ",";
1013                     }
1014                     
1015                     outputString += dupNames[dupNames.size()-1];
1016                     if (tree[node].getBranchLength() != -1) {
1017                         outputString += ":" + toString(0.0);
1018                     }
1019                     
1020                     outputString += ")";
1021                     if (tree[node].getBranchLength() != -1) {
1022                         outputString += ":" + toString(tree[node].getBranchLength());
1023                     }
1024                 }
1025             }else { 
1026                 outputString = tree[node].getName();
1027                 //if there is a branch length then print it
1028                 if (tree[node].getBranchLength() != -1) {
1029                     outputString += ":" + toString(tree[node].getBranchLength());
1030                 }
1031                 
1032                 m->mothurOut("[ERROR]: " + tree[node].getName() + " is not in your namefile, please correct."); m->mothurOutEndLine(); 
1033             }
1034                 
1035             out << outputString;
1036                 }
1037                 
1038         }
1039         catch(exception& e) {
1040                 m->errorOut(e, "Tree", "printBranch");
1041                 exit(1);
1042         }
1043 }
1044 /*****************************************************************/
1045 void Tree::printBranch(int node, ostream& out, string mode) {
1046     try {
1047         
1048         // you are not a leaf
1049                 if (tree[node].getLChild() != -1) {
1050                         out << "(";
1051                         printBranch(tree[node].getLChild(), out, mode);
1052                         out << ",";
1053                         printBranch(tree[node].getRChild(), out, mode);
1054                         out << ")";
1055                         if (mode == "branch") {
1056                                 //if there is a branch length then print it
1057                                 if (tree[node].getBranchLength() != -1) {
1058                                         out << ":" << tree[node].getBranchLength();
1059                                 }
1060                         }else if (mode == "boot") {
1061                                 //if there is a label then print it
1062                                 if (tree[node].getLabel() != -1) {
1063                                         out << tree[node].getLabel();
1064                                 }
1065                         }else if (mode == "both") {
1066                                 if (tree[node].getLabel() != -1) {
1067                                         out << tree[node].getLabel();
1068                                 }
1069                                 //if there is a branch length then print it
1070                                 if (tree[node].getBranchLength() != -1) {
1071                                         out << ":" << tree[node].getBranchLength();
1072                                 }
1073                         }
1074                 }else { //you are a leaf
1075                         string leafGroup = tmap->getGroup(tree[node].getName());
1076                         
1077                         if (mode == "branch") {
1078                                 out << leafGroup; 
1079                                 //if there is a branch length then print it
1080                                 if (tree[node].getBranchLength() != -1) {
1081                                         out << ":" << tree[node].getBranchLength();
1082                                 }
1083                         }else if (mode == "boot") {
1084                                 out << leafGroup; 
1085                                 //if there is a label then print it
1086                                 if (tree[node].getLabel() != -1) {
1087                                         out << tree[node].getLabel();
1088                                 }
1089                         }else if (mode == "both") {
1090                                 out << tree[node].getName();
1091                                 if (tree[node].getLabel() != -1) {
1092                                         out << tree[node].getLabel();
1093                                 }
1094                                 //if there is a branch length then print it
1095                                 if (tree[node].getBranchLength() != -1) {
1096                                         out << ":" << tree[node].getBranchLength();
1097                                 }
1098                         }
1099                 }
1100                 
1101         }
1102         catch(exception& e) {
1103                 m->errorOut(e, "Tree", "printBranch");
1104                 exit(1);
1105         }
1106 }
1107 /*****************************************************************/
1108 void Tree::printBranch(int node, ostream& out, string mode, vector<Node>& theseNodes) {
1109         try {
1110                 
1111                 // you are not a leaf
1112                 if (theseNodes[node].getLChild() != -1) {
1113                         out << "(";
1114                         printBranch(theseNodes[node].getLChild(), out, mode);
1115                         out << ",";
1116                         printBranch(theseNodes[node].getRChild(), out, mode);
1117                         out << ")";
1118                         if (mode == "branch") {
1119                                 //if there is a branch length then print it
1120                                 if (theseNodes[node].getBranchLength() != -1) {
1121                                         out << ":" << theseNodes[node].getBranchLength();
1122                                 }
1123                         }else if (mode == "boot") {
1124                                 //if there is a label then print it
1125                                 if (theseNodes[node].getLabel() != -1) {
1126                                         out << theseNodes[node].getLabel();
1127                                 }
1128                         }else if (mode == "both") {
1129                                 if (theseNodes[node].getLabel() != -1) {
1130                                         out << theseNodes[node].getLabel();
1131                                 }
1132                                 //if there is a branch length then print it
1133                                 if (theseNodes[node].getBranchLength() != -1) {
1134                                         out << ":" << theseNodes[node].getBranchLength();
1135                                 }
1136                         }
1137                 }else { //you are a leaf
1138                         string leafGroup = tmap->getGroup(theseNodes[node].getName());
1139                         
1140                         if (mode == "branch") {
1141                                 out << leafGroup; 
1142                                 //if there is a branch length then print it
1143                                 if (theseNodes[node].getBranchLength() != -1) {
1144                                         out << ":" << theseNodes[node].getBranchLength();
1145                                 }
1146                         }else if (mode == "boot") {
1147                                 out << leafGroup; 
1148                                 //if there is a label then print it
1149                                 if (theseNodes[node].getLabel() != -1) {
1150                                         out << theseNodes[node].getLabel();
1151                                 }
1152                         }else if (mode == "both") {
1153                                 out << theseNodes[node].getName();
1154                                 if (theseNodes[node].getLabel() != -1) {
1155                                         out << theseNodes[node].getLabel();
1156                                 }
1157                                 //if there is a branch length then print it
1158                                 if (theseNodes[node].getBranchLength() != -1) {
1159                                         out << ":" << theseNodes[node].getBranchLength();
1160                                 }
1161                         }
1162                 }
1163                 
1164         }
1165         catch(exception& e) {
1166                 m->errorOut(e, "Tree", "printBranch");
1167                 exit(1);
1168         }
1169 }
1170 /*****************************************************************/
1171
1172 void Tree::printTree() {
1173         
1174         for(int i=0;i<numNodes;i++){
1175                 cout << i << '\t';
1176                 tree[i].printNode();
1177         }
1178         
1179 }
1180
1181 /*****************************************************************/
1182 //this code is a mess and should be rethought...-slw
1183 void Tree::parseTreeFile() {
1184         
1185         //only takes names from the first tree and assumes that all trees use the same names.
1186         try {
1187                 string filename = m->getTreeFile();
1188                 ifstream filehandle;
1189                 m->openInputFile(filename, filehandle);
1190                 int c, comment;
1191                 comment = 0;
1192                 int done = 1;
1193                 
1194                 //ifyou are not a nexus file 
1195                 if((c = filehandle.peek()) != '#') {  
1196                         while((c = filehandle.peek()) != ';') { 
1197                                 while ((c = filehandle.peek()) != ';') {
1198                                         // get past comments
1199                                         if(c == '[') {
1200                                                 comment = 1;
1201                                         }
1202                                         if(c == ']'){
1203                                                 comment = 0;
1204                                         }
1205                                         if((c == '(') && (comment != 1)){ break; }
1206                                         filehandle.get();
1207                                 }
1208
1209                                 done = readTreeString(filehandle); 
1210                                 if (done == 0) { break; }
1211                         }
1212                 //ifyou are a nexus file
1213                 }else if((c = filehandle.peek()) == '#') {
1214                         string holder = "";
1215                                         
1216                         // get past comments
1217                         while(holder != "translate" && holder != "Translate"){  
1218                                 if(holder == "[" || holder == "[!"){
1219                                         comment = 1;
1220                                 }
1221                                 if(holder == "]"){
1222                                         comment = 0;
1223                                 }
1224                                 filehandle >> holder; 
1225
1226                                 //if there is no translate then you must read tree string otherwise use translate to get names
1227                                 if((holder == "tree") && (comment != 1)){       
1228                                         //pass over the "tree rep.6878900 = "
1229                                         while (((c = filehandle.get()) != '(') && ((c = filehandle.peek()) != EOF)) {;}
1230
1231                                         if(c == EOF) { break; }
1232                                         filehandle.putback(c);  //put back first ( of tree.
1233                                         done = readTreeString(filehandle);
1234         
1235                                         break;
1236                                 }
1237                         
1238                                 if (done == 0) { break;  }
1239                         }
1240                         
1241                         //use nexus translation rather than parsing tree to save time
1242                         if((holder == "translate") || (holder == "Translate")) {
1243
1244                                 string number, name, h;
1245                                 h = ""; // so it enters the loop the first time
1246                                 while((h != ";") && (number != ";")) { 
1247                                         filehandle >> number;
1248                                         filehandle >> name;
1249         
1250                                         //c = , until done with translation then c = ;
1251                                         h = name.substr(name.length()-1, name.length()); 
1252                                         name.erase(name.end()-1);  //erase the comma
1253                                         m->Treenames.push_back(number);
1254                                 }
1255                                 if(number == ";") { m->Treenames.pop_back(); }  //in case ';' from translation is on next line instead of next to last name
1256                         }
1257                 }
1258                 filehandle.close();
1259                 
1260                 //for (int i = 0; i < globaldata->Treenames.size(); i++) {
1261 //cout << globaldata->Treenames[i] << endl; }
1262 //cout << globaldata->Treenames.size() << endl;
1263         }
1264         catch(exception& e) {
1265                 m->errorOut(e, "Tree", "parseTreeFile");
1266                 exit(1);
1267         }
1268 }
1269 /*******************************************************/
1270
1271 /*******************************************************/
1272 int Tree::readTreeString(ifstream& filehandle)  {
1273         try {
1274                 int c;
1275                 string name;  //, k
1276                 
1277                 while((c = filehandle.peek()) != ';') { 
1278 //k = c;
1279 //cout << " at beginning of while " <<  k << endl;                      
1280                         if(c == ')')  {    
1281                                 //to pass over labels in trees
1282                                 c=filehandle.get();
1283                                 while((c!=',') && (c != -1) && (c!= ':') && (c!=';')){ c=filehandle.get(); }
1284                                 filehandle.putback(c);
1285                         }
1286                         if(c == ';') { return 0; }
1287                         if(c == -1) { return 0; }
1288                         //if you are a name
1289                         if((c != '(') && (c != ')') && (c != ',') && (c != ':') && (c != '\n') && (c != '\t') && (c != 32)) { //32 is space
1290                                 name = "";
1291                                 c = filehandle.get();
1292                         //k = c;
1293 //cout << k << endl;
1294                                 while ((c != '(') && (c != ')') && (c != ',') && (c != ':') && (c != '\n') && (c != 32) && (c != '\t')) {                       
1295                                         name += c;
1296                                         c = filehandle.get();
1297                         //k = c;
1298 //cout << " in name while " << k << endl;
1299                                 }
1300                                 
1301 //cout << "name = " << name << endl;
1302                                 m->Treenames.push_back(name);
1303                                 filehandle.putback(c);
1304 //k = c;
1305 //cout << " after putback" <<  k << endl;
1306                         } 
1307                         
1308                         if(c  == ':') { //read until you reach the end of the branch length
1309                                 while ((c != '(') && (c != ')') && (c != ',') && (c != ';') && (c != '\n') && (c != '\t') && (c != 32)) {
1310                                         c = filehandle.get();
1311         //k = c;
1312         //cout << " in branch while " << k << endl;
1313                                 }
1314                                 filehandle.putback(c);
1315                         }
1316                 
1317                         c = filehandle.get();
1318 //k = c;
1319         //cout << " here after get " << k << endl;
1320                         if(c == ';') { return 0; }
1321                         if(c == ')') { filehandle.putback(c); }
1322         //k = c;
1323 //cout << k << endl;
1324
1325                 }
1326                 return 0;
1327         }
1328         catch(exception& e) {
1329                 m->errorOut(e, "Tree", "readTreeString");
1330                 exit(1);
1331         }
1332 }       
1333
1334 /*******************************************************/
1335
1336 /*******************************************************/
1337