fixed parsimony for groups

diff --git a/parsimony.cpp b/parsimony.cpp
index 3ef141ed8518eb4ff9a0057152484d2ef85b3968..54960cc516047f0f43ff8cde28f697512d9d6c91 100644 (file)
--- a/parsimony.cpp
+++ b/parsimony.cpp
@@ -17,7 +17,12 @@ EstOutput Parsimony::getValues(Tree* t) {
                data.resize(1,0);
                        
                int score = 0;
-                                       
+               
+               //create pgroups that reflect the groups the user want to use
+               for(int i=t->getNumLeaves();i<t->getNumNodes();i++){
+                       t->tree[i].pGroups = (t->mergeUserGroups(i));
+               }
+
                for(int i=t->getNumLeaves();i<t->getNumNodes();i++){
                        int lc = t->tree[i].getLChild();
                        int rc = t->tree[i].getRChild();
@@ -25,21 +30,25 @@ EstOutput Parsimony::getValues(Tree* t) {
                        int iSize = t->tree[i].pGroups.size();
                        int rcSize = t->tree[rc].pGroups.size();
                        int lcSize = t->tree[lc].pGroups.size();
-                       
+cout << " i groups ";                  
                        //add in all the groups the users wanted
                        for (it = t->tree[i].pGroups.begin(); it != t->tree[i].pGroups.end(); it++) {
-                               if (inUsersGroups(it->first, globaldata->Groups) != true) {  iSize--;  }
+cout << it->first << " ";
+                       //      if (inUsersGroups(it->first, globaldata->Groups) != true) {  iSize--;  }
                        }
+cout << endl << " rc groups ";
                        //add in all the groups the users wanted
                        for (it = t->tree[rc].pGroups.begin(); it != t->tree[rc].pGroups.end(); it++) {
-                               if (inUsersGroups(it->first, globaldata->Groups) != true) {  rcSize--;  }
+cout << it->first << " ";
+                               //if (inUsersGroups(it->first, globaldata->Groups) != true) {  rcSize--;  }
                        }
-                       
+cout << endl << " lc groups ";                 
                        //add in all the groups the users wanted
                        for (it = t->tree[lc].pGroups.begin(); it != t->tree[lc].pGroups.end(); it++) {
-                               if (inUsersGroups(it->first, globaldata->Groups) != true) {  lcSize--;  }
+cout << it->first << " ";
+                               //if (inUsersGroups(it->first, globaldata->Groups) != true) {  lcSize--;  }
                        }
-                       
+cout << endl;                  
                        //if isize are 0 then that branch is to be ignored
                        if (iSize == 0) { }
                        else if ((rcSize == 0) || (lcSize == 0)) { }

diff --git a/parsimonycommand.cpp b/parsimonycommand.cpp
index 554d24bd93f76c94a223d33ce24159babd0f3e78..59dca88751edfb115cc3c83fcc5c85935b52b9d3 100644 (file)
--- a/parsimonycommand.cpp
+++ b/parsimonycommand.cpp
@@ -62,10 +62,13 @@ int ParsimonyCommand::execute() {
                outDist << "RandomTree#" << '\t' << "ParsScore" << endl;
                
                if (randomtree == "") {
+                       copyUserTree = new Tree();
                        //get pscores for users trees
                        for (int i = 0; i < T.size(); i++) {
+                               //copy users tree so that you can redo pgroups 
+                               copyUserTree->getCopy(T[i]);
                                cout << "Processing tree " << i+1 << endl;
-                               userData = pars->getValues(T[i]);  //userData[0] = pscore
+                               userData = pars->getValues(copyUserTree);  //userData[0] = pscore
                                cout << "Tree " << i+1 << " parsimony score = " << userData[0] << endl;
                                //update uscoreFreq
                                it = uscoreFreq.find(userData[0]);

diff --git a/parsimonycommand.h b/parsimonycommand.h
index 326c0c9170ceae07e7061d31c46a25653f52cb7d..2881dc61e200540588691de3f205420139510561 100644 (file)
--- a/parsimonycommand.h
+++ b/parsimonycommand.h
@@ -28,6 +28,7 @@ class ParsimonyCommand : public Command {
                GlobalData* globaldata;
                vector<Tree*> T;           //user trees
                Tree* randT;  //random tree
+               Tree* copyUserTree; 
                TreeMap* tmap; 
                TreeMap* savetmap;
                Parsimony* pars;

diff --git a/tree.cpp b/tree.cpp
index 2ad46f044f71556a9726c93a96d4fcd5ae42a6b6..d332f4771ba0b3f890b407b8f5a0ca31834b46ad 100644 (file)
--- a/tree.cpp
+++ b/tree.cpp
@@ -154,11 +154,79 @@ void Tree::getCopy(Tree* copy) {
 //and you get p[white] = 2, p[black] = 1, but you erase the p[black] because you have a p value higher than 1.
 
 map<string, int> Tree::mergeGroups(int i) {
+       try {
+               int lc = tree[i].getLChild();
+               int rc = tree[i].getRChild();
+               
+               //set parsimony groups to left child
+               map<string,int> parsimony = tree[lc].pGroups;
+               
+               int maxPars = 1;
+
+               //look at right child groups and update maxPars if right child has something higher for that group.
+               for(it=tree[rc].pGroups.begin();it!=tree[rc].pGroups.end();it++){
+                       it2 = parsimony.find(it->first);
+                       if (it2 != parsimony.end()) {
+                               parsimony[it->first]++;
+                       }else {
+                               parsimony[it->first] = 1;
+                       }
+                       
+                       if(parsimony[it->first] > maxPars){
+                               maxPars = parsimony[it->first];
+                       }
+               }
+       
+               // this is true if right child had a greater parsimony for a certain group
+               if(maxPars > 1){
+                       //erase all the groups that are only 1 because you found something with 2.
+                       for(it=parsimony.begin();it!=parsimony.end();it++){
+                               if(it->second == 1){
+                                       parsimony.erase(it->first);
+                                       it--;
+                               }
+                       }
+                       //set one remaining groups to 1
+                       //so with our above example p[white] = 2 would be left and it would become p[white] = 1
+                       for(it=parsimony.begin();it!=parsimony.end();it++){
+                               parsimony[it->first] = 1;
+                       }
+               
+               }
+       
+               return parsimony;
+       }
+       catch(exception& e) {
+               cout << "Standard Error: " << e.what() << " has occurred in the Tree class Function mergeGroups. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
+               exit(1);
+       }
+       catch(...) {
+               cout << "An unknown error has occurred in the Tree class function mergeGroups. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
+               exit(1);
+       }               
+}
+/*****************************************************************/
+//returns a map with a groupname and the number of times that group was seen in the children
+//for instance if your children are white and black then it would return a map with 2 entries
+// p[white] = 1 and p[black] = 1.  Now go up a level and merge that with a node who has p[white] = 1
+//and you get p[white] = 2, p[black] = 1, but you erase the p[black] because you have a p value higher than 1.
+
+map<string, int> Tree::mergeUserGroups(int i) {
        try {
        
                int lc = tree[i].getLChild();
                int rc = tree[i].getRChild();
                
+               //loop through nodes groups removing the ones the user doesn't want
+               for (it = tree[lc].pGroups.begin(); it != tree[lc].pGroups.end(); it++) {
+                       if (inUsersGroups(it->first, globaldata->Groups) != true) { tree[lc].pGroups.erase(it->first); }
+               }
+               
+               //loop through nodes groups removing the ones the user doesn't want
+               for (it = tree[rc].pGroups.begin(); it != tree[rc].pGroups.end(); it++) {
+                       if (inUsersGroups(it->first, globaldata->Groups) != true) { tree[rc].pGroups.erase(it->first); }
+               }
+
                //set parsimony groups to left child
                map<string,int> parsimony = tree[lc].pGroups;
                
@@ -166,7 +234,12 @@ map<string, int> Tree::mergeGroups(int i) {
 
                //look at right child groups and update maxPars if right child has something higher for that group.
                for(it=tree[rc].pGroups.begin();it!=tree[rc].pGroups.end();it++){
-                       parsimony[it->first]++;
+                       it2 = parsimony.find(it->first);
+                       if (it2 != parsimony.end()) {
+                               parsimony[it->first]++;
+                       }else {
+                               parsimony[it->first] = 1;
+                       }
                        
                        if(parsimony[it->first] > maxPars){
                                maxPars = parsimony[it->first];
@@ -202,6 +275,7 @@ map<string, int> Tree::mergeGroups(int i) {
        }               
 }
 
+
 /**************************************************************************************************/
 
 map<string,int> Tree::mergeGcounts(int position) {

diff --git a/tree.h b/tree.h
index 7157a6f0e43993e8278591d3422caee1eeaa8728..55d3c92aabca1083fd300896a813f41374eca860 100644 (file)
--- a/tree.h
+++ b/tree.h
@@ -31,6 +31,7 @@ class Tree {
                void setIndex(string, int);
                int getNumNodes() { return numNodes; }
                int getNumLeaves(){     return numLeaves; }
+               map<string, int> mergeUserGroups(int);  //returns a map with a groupname and the number of times that group was seen in the children
                
                //this function takes the leaf info and populates the non leaf nodes
                void assembleTree();            
@@ -43,8 +44,9 @@ class Tree {
                ofstream out;
                string filename;
                
-               map<string, int>::iterator it;
+               map<string, int>::iterator it, it2;
                map<string, int> mergeGroups(int);  //returns a map with a groupname and the number of times that group was seen in the children
+               
                map<string,int> Tree::mergeGcounts(int);
                void randomTopology();
                void randomBlengths();