pat's fixes to pca command

[mothur.git] / linearalgebra.cpp

diff --git a/linearalgebra.cpp b/linearalgebra.cpp
index 9308627d1335920ea49f9b4c5bbd543de629ae85..e19861912d9c7fc78f9b2cbbb87ecd76e9388ba9 100644 (file)
--- 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;
+       
+}
 
+/*********************************************************************************************************************************/