X-Git-Url: https://git.donarmstrong.com/?a=blobdiff_plain;f=linearalgebra.cpp;h=e19861912d9c7fc78f9b2cbbb87ecd76e9388ba9;hb=f0a594f6676ef5a52d1f122b6de70de2fda08c81;hp=9308627d1335920ea49f9b4c5bbd543de629ae85;hpb=5553e33be3a45eee6bed2ac9a5c4ca0aa0e8d5e4;p=mothur.git

diff --git a/linearalgebra.cpp b/linearalgebra.cpp
index 9308627..e198619 100644
--- a/linearalgebra.cpp
+++ b/linearalgebra.cpp
@@ -53,6 +53,39 @@ vector<vector<double> > LinearAlgebra::matrix_mult(vector<vector<double> > first
 
 /*********************************************************************************************************************************/
 
+void LinearAlgebra::recenter(double offset, vector<vector<double> > D, vector<vector<double> >& G){
+	try {
+		int rank = D.size();
+		
+		vector<vector<double> > A(rank);
+		vector<vector<double> > C(rank);
+		for(int i=0;i<rank;i++){
+			A[i].resize(rank);
+			C[i].resize(rank);
+		}
+		
+		double scale = -1.0000 / (double) rank;
+		
+		for(int i=0;i<rank;i++){
+			A[i][i] = 0.0000;
+			C[i][i] = 1.0000 + scale;
+			for(int j=i+1;j<rank;j++){
+				A[i][j] = A[j][i] = -0.5 * D[i][j] * D[i][j] + offset;
+				C[i][j] = C[j][i] = scale;
+			}
+		}
+		
+		A = matrix_mult(C,A);
+		G = matrix_mult(A,C);
+	}
+	catch(exception& e) {
+		m->errorOut(e, "LinearAlgebra", "recenter");
+		exit(1);
+	}
+	
+}
+/*********************************************************************************************************************************/
+
 //  This function is taken from Numerical Recipes in C++ by Press et al., 2nd edition, pg. 479
 
 int LinearAlgebra::tred2(vector<vector<double> >& a, vector<double>& d, vector<double>& e){
@@ -637,7 +670,7 @@ double LinearAlgebra::calcKendall(vector< vector<double> >& euclidDists, vector<
 			int numWithLowerRank = 0;
 			float thisrank = user[l].score;
 					
-			for (int u = l; u < scores.size(); u++) {
+			for (int u = l+1; u < scores.size(); u++) {
 				if (user[u].score > thisrank) { numWithHigherRank++; }
 				else if (user[u].score < thisrank) { numWithLowerRank++; }
 				count++;
@@ -647,9 +680,6 @@ double LinearAlgebra::calcKendall(vector< vector<double> >& euclidDists, vector<
 			numDisCoor += numWithLowerRank;
 		}
 				
-		//comparing to yourself
-		count -= userDists.size();
-				
 		r = (numCoor - numDisCoor) / (float) count;
 		
 		//divide by zero error
@@ -666,4 +696,30 @@ double LinearAlgebra::calcKendall(vector< vector<double> >& euclidDists, vector<
 
 /*********************************************************************************************************************************/
 
+vector<vector<double> > LinearAlgebra::getObservedEuclideanDistance(vector<vector<double> >& relAbundData){
+
+	int numSamples = relAbundData.size();
+	int numOTUs = relAbundData[0].size();
+	
+	vector<vector<double> > dMatrix(numSamples);
+	for(int i=0;i<numSamples;i++){
+		dMatrix[i].resize(numSamples);
+	}
+	
+	for(int i=0;i<numSamples;i++){
+		for(int j=0;j<numSamples;j++){
+			
+			double d = 0;
+			for(int k=0;k<numOTUs;k++){
+				d += pow((relAbundData[i][k] - relAbundData[j][k]), 2.0000);
+			}
+			dMatrix[i][j] = pow(d, 0.50000);
+			dMatrix[j][i] = dMatrix[i][j];
+			
+		}
+	}
+	return dMatrix;
+	
+}
 
+/*********************************************************************************************************************************/