From 010638fc7580c18a51f5ed650e9ee7589c8fb2d4 Mon Sep 17 00:00:00 2001 From: westcott Date: Mon, 23 Feb 2009 20:09:07 +0000 Subject: [PATCH] fixed parsimony for groups --- parsimony.cpp | 23 ++++++++++---- parsimonycommand.cpp | 5 ++- parsimonycommand.h | 1 + tree.cpp | 76 +++++++++++++++++++++++++++++++++++++++++++- tree.h | 4 ++- 5 files changed, 99 insertions(+), 10 deletions(-) diff --git a/parsimony.cpp b/parsimony.cpp index 3ef141e..54960cc 100644 --- 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();igetNumNodes();i++){ + t->tree[i].pGroups = (t->mergeUserGroups(i)); + } + for(int i=t->getNumLeaves();igetNumNodes();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 554d24b..59dca88 100644 --- 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 326c0c9..2881dc6 100644 --- a/parsimonycommand.h +++ b/parsimonycommand.h @@ -28,6 +28,7 @@ class ParsimonyCommand : public Command { GlobalData* globaldata; vector T; //user trees Tree* randT; //random tree + Tree* copyUserTree; TreeMap* tmap; TreeMap* savetmap; Parsimony* pars; diff --git a/tree.cpp b/tree.cpp index 2ad46f0..d332f47 100644 --- 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 Tree::mergeGroups(int i) { + try { + int lc = tree[i].getLChild(); + int rc = tree[i].getRChild(); + + //set parsimony groups to left child + map 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 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 parsimony = tree[lc].pGroups; @@ -166,7 +234,12 @@ map 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 Tree::mergeGroups(int i) { } } + /**************************************************************************************************/ map Tree::mergeGcounts(int position) { diff --git a/tree.h b/tree.h index 7157a6f..55d3c92 100644 --- 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 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::iterator it; + map::iterator it, it2; map mergeGroups(int); //returns a map with a groupname and the number of times that group was seen in the children + map Tree::mergeGcounts(int); void randomTopology(); void randomBlengths(); -- 2.39.2