X-Git-Url: https://git.donarmstrong.com/?a=blobdiff_plain;f=coverage.cpp;h=dea8dff84c4d50676a179f486bd4c805c73c3632;hb=29351fe29c017d2225366418973ed0e000ce5293;hp=6b514fbc2e47a40c8312ecfa206bc4d9e04a1e50;hpb=81ff02ef4064fe4dec986ea3110c1cb8e1285fe6;p=mothur.git

diff --git a/coverage.cpp b/coverage.cpp
index 6b514fb..dea8dff 100644
--- a/coverage.cpp
+++ b/coverage.cpp
@@ -23,7 +23,7 @@ Coverage::Coverage() {
 }
 
 //**********************************************************************************************************************
-void Coverage::getValues(FullMatrix* matrix, vector< vector< vector<float> > >& data, vector<float> dist, string mode) {
+void Coverage::getValues(FullMatrix* matrix, vector< vector< vector<float> > >& data, vector<float> dist) {
 	try {
 		vector<float> min;
 		vector<string> groups;
@@ -41,19 +41,13 @@ void Coverage::getValues(FullMatrix* matrix, vector< vector< vector<float> > >&
 		//get the minimums for each comparision
 		/**************************************/
 		int count = 0;
-		int a = 0;
-		int b = 0;
+		
 		for (int i = 0; i < numGroups; i++) {
 			for (int j = 0; j < numGroups; j++) {
 			
-				//is this "box" one hte user wants analyzed?
+				//is this "box" one the user wants analyzed?
 				if ((inUsersGroups(globaldata->gGroupmap->namesOfGroups[i], globaldata->Groups) == true) && (inUsersGroups(globaldata->gGroupmap->namesOfGroups[j], globaldata->Groups) == true)) {
 					
-					if (mode == "random") {
-						//create random matrix for this comparison
-						matrix->shuffle(globaldata->gGroupmap->namesOfGroups[i], globaldata->gGroupmap->namesOfGroups[j]);
-					}
-			
 					min = matrix->getMins(count);  //returns vector of mins for "box" requested ie. groups A, B, 0 = AA, 1 = AB, 2 = BA, 3 = BB;
 
 					//find the coverage at this distance
@@ -72,25 +66,109 @@ void Coverage::getValues(FullMatrix* matrix, vector< vector< vector<float> > >&
 						if (index == -1) { index = min.size();  }
 					
 						//save value in data
-						data[k][a][b] = 1.0 - ((min.size()-index)/(float)min.size());
+						data[k][i][j] = 1.0 - ((min.size()-index)/(float)min.size());
 	
 					}
+				}
+				count++;
+			}
+		}
+		
+	}
+	catch(exception& e) {
+		cout << "Standard Error: " << e.what() << " has occurred in the Coverage class Function getValues. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
+		exit(1);
+	}
+	catch(...) {
+		cout << "An unknown error has occurred in the Coverage class function getValues. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
+		exit(1);
+	}	
+}
+
+//**********************************************************************************************************************
+void Coverage::getValues(FullMatrix* matrix, vector< vector< vector<float> > >& data, vector<float> dist, string mode) {
+	try {
+		vector<float> min1;
+		vector<float> min2;
+		vector<string> groups;
+		
+		//initialize data
+		data.resize(dist.size());
+		for (int l = 0; l < data.size(); l++) {
+			data[l].resize(numGroups);
+			for (int k = 0; k < data[l].size(); k++) {
+				data[l][k].push_back(0.0);
+			}
+		}
+		
+		int count = 0;
+		int count2 = 0;
+		
+		//for each box
+		for (int i = 0; i < numGroups; i++) {
+			for (int j = 0; j < numGroups; j++) {
+				
+				if (i != j) {
+					//is this "box" one the user wants analyzed?
+					if ((inUsersGroups(globaldata->gGroupmap->namesOfGroups[i], globaldata->Groups) == true) && (inUsersGroups(globaldata->gGroupmap->namesOfGroups[j], globaldata->Groups) == true)) {
 					
-					//move to next box
-					if (b < numUserGroups-1) {  b++;  }
-					else{ //you are moving to a new row of "boxes"
-						b = 0;
-						a++;
-					}
+						matrix->shuffle(globaldata->gGroupmap->namesOfGroups[i], globaldata->gGroupmap->namesOfGroups[j]);
+					
+						min1 = matrix->getMins(count);  //returns vector of mins for "box" requested ie. groups A, B, 0 = AA, 1 = AB, 2 = BA, 3 = BB;
+						min2 = matrix->getMins(count2);  //returns vector of mins for "box" requested ie. groups A, B, 0 = AA, 1 = AB, 2 = BA, 3 = BB;
+
+						//find the coverage at this distance
+						sort(min1.begin(), min1.end());
+					
+						//find the coverage at this distance
+						sort(min2.begin(), min2.end());
+					
+						float distDiff = 0;
+						
+						//loop through each distance and fill data
+						for (int k = 0; k < data.size(); k++) {
+							//****** coverage of AA **********//
+							int index = -1;
+							//find index in min where value is higher than d
+							for (int m = 0; m < min1.size(); m++) {
+									if (min1[m] > dist[k])	{ index = m; break;	}
+							}
+					
+							// if you don't find one than all the mins are less than d
+							if (index == -1) { index = min1.size();  }
+							
+							//****** coverage of AB **********//
+							int index2 = -1;
+							//find index in min where value is higher than d
+							for (int m = 0; m < min2.size(); m++) {
+									if (min2[m] > dist[k])	{ index2 = m; break;	}
+							}
+					
+							// if you don't find one than all the mins are less than d
+							if (index2 == -1) { index2 = min2.size();  }
 
-					count++;
+							//coverage of ii
+							float covII = 1.0 - ((min1.size()-index)/(float)min1.size());
+							
+							//coverage of ij
+							float covIJ = 1.0 - ((min2.size()-index2)/(float)min2.size());
+							
+							//save value in data (Caa - Cab)^2 * distDiff
+							data[k][i][j] = ((covII-covIJ) * (covII-covIJ)) * distDiff;
+							
+							//update distDiff
+							if (k < data.size() - 1) {
+								distDiff = dist[k+1] - dist[k];	
+							}
+						}
 					
-					if (mode == "random") {
-						//restore matrix to original form for next shuffle
+						//put matrix back to original
 						matrix->restore();
 					}
 				}
+				count2++;
 			}
+			count += numGroups+1; //go from AA to BB to CC
 		}
 		
 	}
@@ -102,7 +180,6 @@ void Coverage::getValues(FullMatrix* matrix, vector< vector< vector<float> > >&
 		cout << "An unknown error has occurred in the Coverage class function getValues. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
 		exit(1);
 	}	
+
 }
-//**********************************************************************************************************************
 
-	
\ No newline at end of file