From: pschloss <pschloss>
Date: Thu, 26 Feb 2009 19:37:22 +0000 (+0000)
Subject: fix weighted unifrac bug in findIndex and randomize labels
X-Git-Url: https://git.donarmstrong.com/?p=mothur.git;a=commitdiff_plain;h=b2c44d69197045c26c1f4e62d248ab882c79a393

fix weighted unifrac bug in findIndex and randomize labels
---

diff --git a/Mothur.xcodeproj/project.pbxproj b/Mothur.xcodeproj/project.pbxproj
index 62d929a..483badd 100644
--- a/Mothur.xcodeproj/project.pbxproj
+++ b/Mothur.xcodeproj/project.pbxproj
@@ -740,6 +740,7 @@
 				GCC_OPTIMIZATION_LEVEL = 3;
 				GCC_WARN_ABOUT_RETURN_TYPE = YES;
 				GCC_WARN_UNUSED_VARIABLE = YES;
+				INSTALL_MODE_FLAG = "a-w,a+rX";
 				PREBINDING = NO;
 				SDKROOT = "$(DEVELOPER_SDK_DIR)/MacOSX10.5.sdk";
 			};
@@ -755,6 +756,7 @@
 				GCC_OPTIMIZATION_LEVEL = 0;
 				GCC_WARN_ABOUT_RETURN_TYPE = YES;
 				GCC_WARN_UNUSED_VARIABLE = YES;
+				INSTALL_MODE_FLAG = "a-w,a+rX";
 				PREBINDING = NO;
 				SDKROOT = "$(DEVELOPER_SDK_DIR)/MacOSX10.5.sdk";
 			};
diff --git a/parsimony.cpp b/parsimony.cpp
index da8412c..3840514 100644
--- a/parsimony.cpp
+++ b/parsimony.cpp
@@ -10,6 +10,7 @@
 #include "parsimony.h"
 
 /**************************************************************************************************/
+
 EstOutput Parsimony::getValues(Tree* t) {
 	try {
 		globaldata = GlobalData::getInstance();
@@ -38,12 +39,8 @@ EstOutput Parsimony::getValues(Tree* t) {
 			else if(iSize > rcSize || iSize > lcSize){
 				score++;
 			}
-			cout << i << ' ' << score << ": ";
-			t->tree[i].printNode();
 		} 
 		
-		//string hold;
-		//cin >> hold;
 		
 		data[0] = score;
 		
diff --git a/tree.cpp b/tree.cpp
index 2689b03..3b3b684 100644
--- a/tree.cpp
+++ b/tree.cpp
@@ -354,33 +354,27 @@ void Tree::randomLabels() {
 
 void Tree::randomLabels(string groupA, string groupB) {
 	try {
-		for(int i = 0; i < numLeaves; i++) {
-			int z;
-			//get random index to switch with
-			z = int((float)(i+1) * (float)(rand()) / ((float)RAND_MAX+1.0));	
-			
-			//you only want to randomize the nodes that are from a group the user wants analyzed, so
-			//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.
-			if (((tree[z].getGroup() == groupA) || (tree[z].getGroup() == groupB)) && ((tree[i].getGroup() == groupA) || (tree[i].getGroup() == groupB))) {
-				//switches node i and node z's info.
-				map<string,int> lib_hold = tree[z].pGroups;
-				tree[z].pGroups = (tree[i].pGroups);
-				tree[i].pGroups = (lib_hold);
-				
-				string zgroup = tree[z].getGroup();
-				tree[z].setGroup(tree[i].getGroup());
-				tree[i].setGroup(zgroup);
-				
-				string zname = tree[z].getName();
-				tree[z].setName(tree[i].getName());
-				tree[i].setName(zname);
+		int numSeqsA = globaldata->gTreemap->seqsPerGroup[groupA];
+		int numSeqsB = globaldata->gTreemap->seqsPerGroup[groupB];
+
+		vector<string> randomGroups(numSeqsA+numSeqsB, groupA);
+		for(int i=numSeqsA;i<randomGroups.size();i++){
+			randomGroups[i] = groupB;
+		}
+		random_shuffle(randomGroups.begin(), randomGroups.end());
 				
-				map<string,int> gcount_hold = tree[z].pcount;
-				tree[z].pcount = (tree[i].pcount);
-				tree[i].pcount = (gcount_hold);
+		int randomCounter = 0;				
+		for(int i=0;i<numLeaves;i++){
+			if(tree[i].getGroup() == groupA || tree[i].getGroup() == groupB){
+				tree[i].setGroup(randomGroups[randomCounter]);
+				tree[i].pcount.clear();
+				tree[i].pcount[randomGroups[randomCounter]] = 1;
+				tree[i].pGroups.clear();
+				tree[i].pGroups[randomGroups[randomCounter]] = 1;
+				randomCounter++;
 			}
 		}
-	}
+	}		
 	catch(exception& e) {
 		cout << "Standard Error: " << e.what() << " has occurred in the Tree class Function randomLabels. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
 		exit(1);
@@ -528,7 +522,7 @@ void Tree::printBranch(int node) {
 			printBranch(tree[node].getRChild());
 			out << ")";
 		}else { //you are a leaf
-			out << tree[node].getName() << ":" << tree[node].getBranchLength();
+			out << tree[node].getGroup() << ":" << tree[node].getBranchLength();
 		}
 		
 	}
diff --git a/unifracweightedcommand.cpp b/unifracweightedcommand.cpp
index 8fc5b3b..a99459d 100644
--- a/unifracweightedcommand.cpp
+++ b/unifracweightedcommand.cpp
@@ -72,10 +72,10 @@ int UnifracWeightedCommand::execute() {
 			
 			//get scores for random trees
 			for (int j = 0; j < iters; j++) {
-				int n = 1;
+//				int n = 1;
 				int count = 0;
 				for (int r=0; r<numGroups; r++) { 
-					for (int l = n; l < numGroups; l++) {
+					for (int l = r+1; l < numGroups; l++) {
 						//copy T[i]'s info.
 						randT->getCopy(T[i]);
 						 
@@ -90,18 +90,18 @@ int UnifracWeightedCommand::execute() {
 							//get wscore of random tree
 							randomData = weighted->getValues(randT, tmap->namesOfGroups[r], tmap->namesOfGroups[l]);
 						}
+//						randT->createNewickFile("hold"+toString(r)+toString(l)+toString(j));
 
 						//save scores
 						rScores[count].push_back(randomData[0]);
 						validScores[count][randomData[0]] = randomData[0];
 						count++;
 					}
-					n++;
+//					n++;
 				}
 			}
 
 			removeValidScoresDuplicates(); 
-			
 			//find the signifigance of the score for summary file
 			for (int f = 0; f < numComp; f++) {
 				//sort random scores
@@ -111,7 +111,7 @@ int UnifracWeightedCommand::execute() {
 				//so if you have 1000 random trees the index returned is 100 
 				//then there are 900 trees with a score greater then you. 
 				//giving you a signifigance of 0.900
-				int index = findIndex(userData[f]);    if (index == -1) { cout << "error in UnifracWeightedCommand" << endl; exit(1); } //error code
+				int index = findIndex(userData[f], f);    if (index == -1) { cout << "error in UnifracWeightedCommand" << endl; exit(1); } //error code
 			
 				//the signifigance is the number of trees with the users score or higher 
 				WScoreSig.push_back((iters-index)/(float)iters);
@@ -228,15 +228,12 @@ void UnifracWeightedCommand::removeValidScoresDuplicates() {
 }
 
 /***********************************************************/
-int UnifracWeightedCommand::findIndex(float score) {
+int UnifracWeightedCommand::findIndex(float score, int index) {
 	try{
-		for (int e = 0; e < numComp; e++) {
-			for (int i = 0; i < rScores[e].size(); i++) {
-//cout << rScores[e][i] << " number " << i << endl;
-				if (rScores[e][i] >= score) { return i; }
-			}
+		for (int i = 0; i < rScores[index].size(); i++) {
+			if (rScores[index][i] >= score)	{	return i;	}
 		}
-		return -1;
+		return rScores[index].size();
 	}
 	catch(exception& e) {
 		cout << "Standard Error: " << e.what() << " has occurred in the UnifracWeightedCommand class Function findIndex. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
diff --git a/unifracweightedcommand.h b/unifracweightedcommand.h
index d309566..e8daaef 100644
--- a/unifracweightedcommand.h
+++ b/unifracweightedcommand.h
@@ -47,7 +47,7 @@ class UnifracWeightedCommand : public Command {
 		void printWSummaryFile();
 	//	void printWeightedFile();  
 		void removeValidScoresDuplicates();
-		int findIndex(float);
+		int findIndex(float, int);
 		void setGroups(); 
 };