From 06da1833191d5dee7e206f436d9631680f1e7c42 Mon Sep 17 00:00:00 2001
From: westcott <westcott>
Date: Mon, 23 Feb 2009 17:48:26 +0000
Subject: [PATCH] fixed parsimony for groups

---
 parsimony.cpp | 36 ++++++++++++------------------------
 readtree.cpp  | 23 +++++++++++++++++++++--
 treemap.h     |  2 ++
 3 files changed, 35 insertions(+), 26 deletions(-)

diff --git a/parsimony.cpp b/parsimony.cpp
index 6239f75..3ef141e 100644
--- a/parsimony.cpp
+++ b/parsimony.cpp
@@ -17,47 +17,35 @@ EstOutput Parsimony::getValues(Tree* t) {
 		data.resize(1,0);
 			
 		int score = 0;
-			
+					
 		for(int i=t->getNumLeaves();i<t->getNumNodes();i++){
 			int lc = t->tree[i].getLChild();
 			int rc = t->tree[i].getRChild();
 			
-			int iSize = 0;
-			int rcSize = 0;
-			int lcSize = 0;
-
+			int iSize = t->tree[i].pGroups.size();
+			int rcSize = t->tree[rc].pGroups.size();
+			int lcSize = t->tree[lc].pGroups.size();
+			
 			//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++;  }
+				if (inUsersGroups(it->first, globaldata->Groups) != true) {  iSize--;  }
 			}
-
-			//if that leaves no groups give it 1 so it will cause no change to parent
-			if (iSize == 0) { iSize++; }
-			
 			//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++;  }
+				if (inUsersGroups(it->first, globaldata->Groups) != true) {  rcSize--;  }
 			}
 			
-			//if that leaves no groups give it 1 so it will cause no change to parent
-			if (rcSize == 0) { rcSize++; }
-
-				
 			//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++;  }
+				if (inUsersGroups(it->first, globaldata->Groups) != true) {  lcSize--;  }
 			}
 			
-			//if that leaves no groups give it 1 so it will cause no change to parent
-			if (lcSize == 0) { lcSize++; }
-
-
+			//if isize are 0 then that branch is to be ignored
+			if (iSize == 0) { }
+			else if ((rcSize == 0) || (lcSize == 0)) { }
 			//if you have more groups than either of your kids then theres been a change.
-			 if(iSize > rcSize || iSize > lcSize){
+			else if(iSize > rcSize || iSize > lcSize){
 				score++;
-
 			}
 		} 
 		
diff --git a/readtree.cpp b/readtree.cpp
index 1910c4a..81a7729 100644
--- a/readtree.cpp
+++ b/readtree.cpp
@@ -331,9 +331,28 @@ int ReadNewickTree::readNewickInt(istream& f, int& n, Tree* T) {
 			//find index in tree of name
 			int n1 = T->getIndex(name);
 			
-			if(n1 == -1){cerr << "Name: " << name << " not found\n"; exit(1);}
+			//adds sequence names that are not in group file to the "xxx" group
+			if(n1 == -1) {
+				cerr << "Name: " << name << " not found in your groupfile and it will be ignored. \n";
+				globaldata->gTreemap->namesOfSeqs.push_back(name);
+				globaldata->gTreemap->treemap[name].groupname = "xxx";
+				globaldata->gTreemap->treemap[name].vectorIndex = (globaldata->gTreemap->namesOfSeqs.size() - 1);
+				
+				map<string, int>::iterator it;
+				it = globaldata->gTreemap->seqsPerGroup.find("xxx");
+				if (it == globaldata->gTreemap->seqsPerGroup.end()) { //its a new group
+					globaldata->gTreemap->namesOfGroups.push_back("xxx");
+					globaldata->gTreemap->seqsPerGroup["xxx"] = 1;
+				}else {
+					globaldata->gTreemap->seqsPerGroup["xxx"]++;
+				}
+				
+				//find index in tree of name
+				n1 = T->getIndex(name);
+				group = "xxx";
+			}
 			
-			else T->tree[n1].setGroup(group);
+			T->tree[n1].setGroup(group);
 		
 			T->tree[n1].setChildren(-1,-1);
 		
diff --git a/treemap.h b/treemap.h
index 0dc74ee..941ca9a 100644
--- a/treemap.h
+++ b/treemap.h
@@ -39,6 +39,7 @@ public:
 	map<string, GroupIndex> treemap; //sequence name and <groupname, vector index>
 	void print(ostream&);
 	
+	
 private:
 	ifstream fileHandle;
 	string groupFileName;
@@ -47,6 +48,7 @@ private:
 	map<string, int>::iterator it2;
 	void setNamesOfGroups(string); 
 	
+	
 };
 
 #endif
-- 
2.39.2