working on pca

diff --git a/linearalgebra.cpp b/linearalgebra.cpp
index 33464ca5b3519f615b1029b22793694d595eee2d..6970c6429481f94cb5d9b7a26bd990407410ff8d 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){

diff --git a/linearalgebra.h b/linearalgebra.h
index e1f6f7ce7ff8e298950a1fb31970405726536a4d..332101853e97ee2ff630f6792dc926c4aa9f4da4 100644 (file)
--- a/linearalgebra.h
+++ b/linearalgebra.h
@@ -20,6 +20,7 @@ public:
        ~LinearAlgebra() {}
        
        vector<vector<double> > matrix_mult(vector<vector<double> >, vector<vector<double> >);
+       void recenter(double, vector<vector<double> >, vector<vector<double> >&);
        int tred2(vector<vector<double> >&, vector<double>&, vector<double>&);
        int qtli(vector<double>&, vector<double>&, vector<vector<double> >&);
        vector< vector<double> > calculateEuclidianDistance(vector<vector<double> >&, int); //pass in axes and number of dimensions

diff --git a/pcacommand.cpp b/pcacommand.cpp
index 3ede67a85e853de3bd5cbd75ad358145e943674d..23c066b7beff9dd0d6c1da089221a07e510411ae 100644 (file)
--- a/pcacommand.cpp
+++ b/pcacommand.cpp
@@ -314,7 +314,7 @@ int PCACommand::process(vector<SharedRAbundFloatVector*>& lookupFloat){
                        out << endl;
                }
                
-               //matrix = linearCalc.matrix_mult(matrix, transposeMatrix);     
+               matrix = linearCalc.matrix_mult(matrix, transposeMatrix);       
                
                out << endl << endl << "matrix mult" << endl;
                for (int i = 0; i < matrix.size(); i++) {

diff --git a/pcacommand.h b/pcacommand.h
index cf4348ca92e393e9e6d46efd7da8f058bab7f422..eedfcf35c58d508cb37ebc41416ca0a3bfc7b1c8 100644 (file)
--- a/pcacommand.h
+++ b/pcacommand.h
@@ -39,7 +39,7 @@ private:
        map<string, vector<string> > outputTypes;
        LinearAlgebra linearCalc;
        
-       vector< vector<double> > createMatrix(vector<SharedRAbundFloatVector*>);
+       //vector< vector<double> > createMatrix(vector<SharedRAbundFloatVector*>);
        int process(vector<SharedRAbundFloatVector*>&);
        void output(string, vector<string>, vector<vector<double> >&, vector<double>);
        

diff --git a/pcoacommand.cpp b/pcoacommand.cpp
index 4be5ffb8c5b63c4202666aecd57ec58c2be05374..3cd844a36d5a5c253786b4eb29745388b042d2d4 100644 (file)
--- a/pcoacommand.cpp
+++ b/pcoacommand.cpp
@@ -169,12 +169,12 @@ int PCOACommand::execute(){
                vector<double> d;
                vector<double> e;
                vector<vector<double> > G = D;
-               vector<vector<double> > copy_G;
+               //vector<vector<double> > copy_G;
                                
                m->mothurOut("\nProcessing...\n\n");
                
                for(int count=0;count<2;count++){
-                       recenter(offset, D, G);                                 if (m->control_pressed) { return 0; }
+                       linearCalc.recenter(offset, D, G);              if (m->control_pressed) { return 0; }
                        linearCalc.tred2(G, d, e);                              if (m->control_pressed) { return 0; }
                        linearCalc.qtli(d, e, G);                               if (m->control_pressed) { return 0; }
                        offset = d[d.size()-1];
@@ -230,40 +230,6 @@ void PCOACommand::get_comment(istream& f, char begin, char end){
 }      
 /*********************************************************************************************************************************/
 
-void PCOACommand::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 = linearCalc.matrix_mult(C,A);
-               G = linearCalc.matrix_mult(A,C);
-       }
-       catch(exception& e) {
-               m->errorOut(e, "PCOACommand", "recenter");
-               exit(1);
-       }
-
-}
-
-/*********************************************************************************************************************************/
-
 void PCOACommand::output(string fnameRoot, vector<string> name_list, vector<vector<double> >& G, vector<double> d) {
        try {
                int rank = name_list.size();

diff --git a/pcoacommand.h b/pcoacommand.h
index 02e562ec23f97d1fdbd0e6e18c949ae3ffc44c54..d923006bdb663bfd8a782a59738681f951eca79f 100644 (file)
--- a/pcoacommand.h
+++ b/pcoacommand.h
@@ -37,7 +37,6 @@ private:
        LinearAlgebra linearCalc;
        
        void get_comment(istream&, char, char);
-       void recenter(double, vector<vector<double> >, vector<vector<double> >&);
        void output(string, vector<string>, vector<vector<double> >&, vector<double>);
        
 };