From 943d3038553f82378567d3b651969146716e5862 Mon Sep 17 00:00:00 2001
From: westcott <westcott>
Date: Tue, 22 Mar 2011 12:43:51 +0000
Subject: [PATCH] working on pca

---
 linearalgebra.cpp | 33 +++++++++++++++++++++++++++++++++
 linearalgebra.h   |  1 +
 pcacommand.cpp    |  2 +-
 pcacommand.h      |  2 +-
 pcoacommand.cpp   | 38 ++------------------------------------
 pcoacommand.h     |  1 -
 6 files changed, 38 insertions(+), 39 deletions(-)

diff --git a/linearalgebra.cpp b/linearalgebra.cpp
index 33464ca..6970c64 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){
diff --git a/linearalgebra.h b/linearalgebra.h
index e1f6f7c..3321018 100644
--- 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 3ede67a..23c066b 100644
--- 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 cf4348c..eedfcf3 100644
--- 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 4be5ffb..3cd844a 100644
--- 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 02e562e..d923006 100644
--- 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>);
 	
 };
-- 
2.39.2