X-Git-Url: https://git.donarmstrong.com/?a=blobdiff_plain;f=homovacommand.cpp;h=a8e2d13f5003f0a8cba0c637e48573aebf134828;hb=5553e33be3a45eee6bed2ac9a5c4ca0aa0e8d5e4;hp=f9bfb1067dd3b1372801c72cf8145b8f8b4e2543;hpb=17d3033f832c6b8fff5875eb4de29f44a926d4dd;p=mothur.git

diff --git a/homovacommand.cpp b/homovacommand.cpp
index f9bfb10..a8e2d13 100644
--- a/homovacommand.cpp
+++ b/homovacommand.cpp
@@ -66,8 +66,7 @@ vector<string> HomovaCommand::getValidParameters(){
 //**********************************************************************************************************************
 HomovaCommand::HomovaCommand(){	
 	try {
-		abort = true;
-		//initialize outputTypes
+		abort = true; calledHelp = true; 
 		vector<string> tempOutNames;
 		outputTypes["homova"] = tempOutNames;
 	}
@@ -105,12 +104,12 @@ vector<string> HomovaCommand::getRequiredFiles(){
 HomovaCommand::HomovaCommand(string option) {
 	try {
 		globaldata = GlobalData::getInstance();
-		abort = false;
+		abort = false; calledHelp = false;   
 		allLines = 1;
 		labels.clear();
 		
 		//allow user to run help
-		if(option == "help") { help(); abort = true; }
+		if(option == "help") { help(); abort = true; calledHelp = true; }
 		
 		else {
 			//valid paramters for this command
@@ -319,6 +318,7 @@ HomovaCommand::HomovaCommand(string option) {
 
 void HomovaCommand::help(){
 	try {
+		m->mothurOut("Referenced: Stewart CN, Excoffier L (1996). Assessing population genetic structure and variability with RAPD data: Application to Vaccinium macrocarpon (American Cranberry). J Evol Biol 9: 153-71.\n");
 		m->mothurOut("The homova command can only be executed after a successful read.otu command of a list and group or shared file, or by providing a phylip formatted distance matrix.\n");
 		m->mothurOut("The homova command outputs a .homova file. \n");
 		m->mothurOut("The homova command parameters are phylip, iters, groups, label, design, sets and processors.  The design parameter is required.\n");
@@ -349,7 +349,7 @@ HomovaCommand::~HomovaCommand(){}
 int HomovaCommand::execute(){
 	try {
 		
-		if (abort == true) { return 0; }
+		if (abort == true) { if (calledHelp) { return 0; }  return 2;	}
 		
 		//read design file
 		designMap = new GroupMap(designfile);
@@ -822,21 +822,17 @@ int HomovaCommand::calcHomova(ofstream& out, int numTreatments, vector<string> t
 		map<string, double> SSWithin;
 		map<string, double>::iterator it;
 		
-		//SSTotal = 0.0;
-		SSTotal = calcTotal(numTreatments);
-		cout << "sstotal = " << SSTotal << " (n-P) = " << (matrix.size() - numTreatments) << endl;
+		SSTotal = calcWithin(matrix, numTreatments, thisCombosLookupSets);
 		int numSamples = matrix.size();
 		
 		//calc BValue
 		map<string, int> counts;
-		SSWithin = calcWithin(matrix, numTreatments, thisCombosLookupSets, counts);
+		SSWithin = calcWithinEach(matrix, numTreatments, thisCombosLookupSets, counts);
 				
 		double sum = 0.0;
 		double sumDenom = 0.0;
 		for (it = SSWithin.begin(); it != SSWithin.end(); it++) {
-			cout << it->first << '\t' << it->second << '\t' << (counts[it->first] - numTreatments) << endl;
 			double temp2 = (it->second) / (double) (counts[it->first] - 1);
-			cout << "sumTerm = " << temp2 << '\t' << "ln sumTerm = " << log(temp2) << endl;
 			sum += ((counts[it->first] - 1) * log(temp2));
 			sumDenom += ((1 / (double) (counts[it->first] - 1)) - ( 1 / (double) (numSamples - numTreatments) ));
 		}
@@ -847,7 +843,6 @@ int HomovaCommand::calcHomova(ofstream& out, int numTreatments, vector<string> t
 		double denom = 1 + ((1 / (double) (3 * (numTreatments - 1))) * sumDenom);
 		
 		BValue = numerator / denom;
-		cout << "numeratorTerm1 = " << numeratorTerm1 << " sum = " << sum << " sumDenom = " << sumDenom << " numerator = " << numerator << " denom = " << denom << " B = " << BValue << endl;	
 		
 		//calc Pvalue
 		int count = 0;
@@ -858,8 +853,10 @@ int HomovaCommand::calcHomova(ofstream& out, int numTreatments, vector<string> t
 			vector<string> copyNames = thisCombosLookupSets;
 			random_shuffle(copyNames.begin(), copyNames.end());
 			
+			SSTotal = calcWithin(matrix, numTreatments, copyNames);
+			
 			counts.clear();
-			map<string, double> randomSSWithin = calcWithin(matrix, numTreatments, copyNames, counts);
+			map<string, double> randomSSWithin = calcWithinEach(matrix, numTreatments, copyNames, counts);
 			
 			sum = 0.0;
 			sumDenom = 0.0;
@@ -869,12 +866,14 @@ int HomovaCommand::calcHomova(ofstream& out, int numTreatments, vector<string> t
 				sumDenom += ((1 / (double) (counts[it->first] - 1)) - ( 1 / (double) (numSamples - numTreatments) ));
 			}
 			
+			temp = SSTotal / (double) (numSamples - numTreatments);
+			numeratorTerm1 = (numSamples - numTreatments) * log(temp);
 			numerator = numeratorTerm1 - sum;
 			denom = 1 + ((1 / (double) (3 * (numTreatments - 1))) * sumDenom);
 			
 			double randomBValue = numerator / denom;
 			
-			if (randomBValue >= BValue) { count++; }
+			if (randomBValue <= BValue) { count++; }
 		}
 				
 		pValue = count / (float) iters;
@@ -889,7 +888,7 @@ int HomovaCommand::calcHomova(ofstream& out, int numTreatments, vector<string> t
 	}
 }
 //**********************************************************************************************************************
-map<string, double> HomovaCommand::calcWithin(vector< vector<double> >& thisMatrix, int numTreatments, vector<string> thisCombosLookupSets, map<string, int>& count) {
+map<string, double> HomovaCommand::calcWithinEach(vector< vector<double> >& thisMatrix, int numTreatments, vector<string> thisCombosLookupSets, map<string, int>& count) {
 	try {
 		map<string, double> within; //maps treatment to within value
 		map<string, double>::iterator it;
@@ -928,29 +927,33 @@ map<string, double> HomovaCommand::calcWithin(vector< vector<double> >& thisMatr
 		return within;
 	}
 	catch(exception& e) {
-		m->errorOut(e, "HomovaCommand", "calcWithin");
+		m->errorOut(e, "HomovaCommand", "calcWithinEach");
 		exit(1);
 	}
 }
 //**********************************************************************************************************************
-double HomovaCommand::calcTotal(int numTreatments) {
+double HomovaCommand::calcWithin(vector< vector<double> >& thisMatrix, int numTreatments, vector<string> thisCombosLookupSets) {
 	try {
 		double total = 0.0;
 		
 		//traverse lower triangle
-		for (int k = 0; k < matrix.size(); k++) { 
-			for (int l = k; l < matrix[k].size(); l++) {
-				total += (matrix[k][l] * matrix[k][l]); //dij^2
+		for (int k = 0; k < thisMatrix.size(); k++) { 
+			for (int l = k; l < thisMatrix[k].size(); l++) {
+				
+				//if you are from the same treatment then eij is 1 so add, else eij = 0
+				if (thisCombosLookupSets[k] == thisCombosLookupSets[l]) { 
+					total += (thisMatrix[k][l] * thisMatrix[k][l]); //dij^2
+				}
 			}
 		}
 		
-		//1 / numSamples
-		total *= (1.0 / (float) matrix.size());
+		//1 / (numSamples / numTreatments)
+		total *= (1.0 / (float) (thisMatrix.size() / (float) numTreatments));
 		
 		return total;
 	}
 	catch(exception& e) {
-		m->errorOut(e, "HomovaCommand", "calcTotal");
+		m->errorOut(e, "HomovaCommand", "calcWithin");
 		exit(1);
 	}
 }