X-Git-Url: https://git.donarmstrong.com/?a=blobdiff_plain;f=pcoacommand.cpp;fp=pcoacommand.cpp;h=9375144bcb342faf364a38eee6eca63226099659;hb=605ab6fa594317a38f0df7bb6797740c735b2348;hp=cb2e758be2ac7aa6ca3afcd5079d2143b59e3ead;hpb=4b6e3f7b5543822a2cca4fb076ab6af2ce8ca62d;p=mothur.git

diff --git a/pcoacommand.cpp b/pcoacommand.cpp
index cb2e758..9375144 100644
--- a/pcoacommand.cpp
+++ b/pcoacommand.cpp
@@ -171,14 +171,13 @@ int PCOACommand::execute(){
 		vector<double> e;
 		vector<vector<double> > G = D;
 		vector<vector<double> > copy_G;
-		//int rank = D.size();
-		
+				
 		m->mothurOut("\nProcessing...\n\n");
 		
 		for(int count=0;count<2;count++){
-			recenter(offset, D, G);		if (m->control_pressed) { return 0; }
-			tred2(G, d, e);				if (m->control_pressed) { return 0; }
-			qtli(d, e, G);				if (m->control_pressed) { return 0; }
+			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];
 			if(offset > 0.0) break;
 		} 
@@ -327,13 +326,6 @@ double PCOACommand::calcPearson(vector< vector<double> >& euclidDists, vector< v
 }
 /*********************************************************************************************************************************/
 
-inline double SIGN(const double a, const double b)
-{
-    return b>=0 ? (a>=0 ? a:-a) : (a>=0 ? -a:a);
-}
-
-/*********************************************************************************************************************************/
-
 void PCOACommand::get_comment(istream& f, char begin, char end){
 	try {
 		char d=f.get();
@@ -444,39 +436,6 @@ void PCOACommand::read(string fname, vector<string>& names, vector<vector<double
 
 /*********************************************************************************************************************************/
 
-double PCOACommand::pythag(double a, double b)	{	return(pow(a*a+b*b,0.5));	}
-
-/*********************************************************************************************************************************/
-
-void PCOACommand::matrix_mult(vector<vector<double> > first, vector<vector<double> > second, vector<vector<double> >& product){
-	try {
-		int first_rows = first.size();
-		int first_cols = first[0].size();
-		int second_cols = second[0].size();
-		
-		product.resize(first_rows);
-		for(int i=0;i<first_rows;i++){
-			product[i].resize(second_cols);
-		}
-		
-		for(int i=0;i<first_rows;i++){
-			for(int j=0;j<second_cols;j++){
-				product[i][j] = 0.0;
-				for(int k=0;k<first_cols;k++){
-					product[i][j] += first[i][k] * second[k][j];
-				}
-			}
-		}
-	}
-	catch(exception& e) {
-		m->errorOut(e, "PCOACommand", "matrix_mult");
-		exit(1);
-	}
-
-}
-
-/*********************************************************************************************************************************/
-
 void PCOACommand::recenter(double offset, vector<vector<double> > D, vector<vector<double> >& G){
 	try {
 		int rank = D.size();
@@ -499,8 +458,8 @@ void PCOACommand::recenter(double offset, vector<vector<double> > D, vector<vect
 			}
 		}
 		
-		matrix_mult(C,A,A);
-		matrix_mult(A,C,G);
+		A = linearCalc.matrix_mult(C,A);
+		G = linearCalc.matrix_mult(A,C);
 	}
 	catch(exception& e) {
 		m->errorOut(e, "PCOACommand", "recenter");
@@ -511,182 +470,6 @@ void PCOACommand::recenter(double offset, vector<vector<double> > D, vector<vect
 
 /*********************************************************************************************************************************/
 
-//  This function is taken from Numerical Recipes in C++ by Press et al., 2nd edition, pg. 479
-
-void PCOACommand::tred2(vector<vector<double> >& a, vector<double>& d, vector<double>& e){
-	try {
-		double scale, hh, h, g, f;
-		
-		int n = a.size();
-		
-		d.resize(n);
-		e.resize(n);
-		
-		for(int i=n-1;i>0;i--){
-			int l=i-1;
-			h = scale = 0.0000;
-			if(l>0){
-				for(int k=0;k<l+1;k++){
-					scale += fabs(a[i][k]);
-				}
-				if(scale == 0.0){
-					e[i] = a[i][l];
-				}
-				else{
-					for(int k=0;k<l+1;k++){
-						a[i][k] /= scale;
-						h += a[i][k] * a[i][k];
-					}
-					f = a[i][l];
-					g = (f >= 0.0 ? -sqrt(h) : sqrt(h));
-					e[i] = scale * g;
-					h -= f * g;
-					a[i][l] = f - g;
-					f = 0.0;
-					for(int j=0;j<l+1;j++){
-						a[j][i] = a[i][j] / h;
-						g = 0.0;
-						for(int k=0;k<j+1;k++){
-							g += a[j][k] * a[i][k];
-						}
-						for(int k=j+1;k<l+1;k++){
-							g += a[k][j] * a[i][k];
-						}
-						e[j] = g / h;
-						f += e[j] * a[i][j];
-					}
-					hh = f / (h + h);
-					for(int j=0;j<l+1;j++){
-						f = a[i][j];
-						e[j] = g = e[j] - hh * f;
-						for(int k=0;k<j+1;k++){
-							a[j][k] -= (f * e[k] + g * a[i][k]);
-						}
-					}
-				}
-			}
-			else{
-				e[i] = a[i][l];
-			}
-			
-			d[i] = h;
-		}
-		
-		d[0] = 0.0000;
-		e[0] = 0.0000;
-		
-		for(int i=0;i<n;i++){
-			int l = i;
-			if(d[i] != 0.0){
-				for(int j=0;j<l;j++){
-					g = 0.0000;
-					for(int k=0;k<l;k++){
-						g += a[i][k] * a[k][j];
-					}
-					for(int k=0;k<l;k++){
-						a[k][j] -= g * a[k][i];
-					}
-				}
-			}
-			d[i] = a[i][i];
-			a[i][i] = 1.0000;
-			for(int j=0;j<l;j++){
-				a[j][i] = a[i][j] = 0.0;
-			}
-		}
-	}
-	catch(exception& e) {
-		m->errorOut(e, "PCOACommand", "tred2");
-		exit(1);
-	}
-
-}
-
-/*********************************************************************************************************************************/
-
-//  This function is taken from Numerical Recipes in C++ by Press et al., 2nd edition, pg. 479
-
-void PCOACommand::qtli(vector<double>& d, vector<double>& e, vector<vector<double> >& z) {
-	try {
-		int m, i, iter;
-		double s, r, p, g, f, dd, c, b;
-		
-		int n = d.size();
-		for(int i=1;i<=n;i++){
-			e[i-1] = e[i];
-		}
-		e[n-1] = 0.0000;
-		
-		for(int l=0;l<n;l++){
-			iter = 0;
-			do {
-				for(m=l;m<n-1;m++){
-					dd = fabs(d[m]) + fabs(d[m+1]);
-					if(fabs(e[m])+dd == dd) break;
-				}
-				if(m != l){
-					if(iter++ == 30) cerr << "Too many iterations in tqli\n";
-					g = (d[l+1]-d[l]) / (2.0 * e[l]);
-					r = pythag(g, 1.0);
-					g = d[m] - d[l] + e[l] / (g + SIGN(r,g));
-					s = c = 1.0;
-					p = 0.0000;
-					for(i=m-1;i>=l;i--){
-						f = s * e[i];
-						b = c * e[i];
-						e[i+1] = (r=pythag(f,g));
-						if(r==0.0){
-							d[i+1] -= p;
-							e[m] = 0.0000;
-							break;
-						}
-						s = f / r;
-						c = g / r;
-						g = d[i+1] - p;
-						r = (d[i] - g) * s + 2.0 * c * b;
-						d[i+1] = g + ( p = s * r);
-						g = c * r - b;
-						for(int k=0;k<n;k++){
-							f = z[k][i+1];
-							z[k][i+1] = s * z[k][i] + c * f;
-							z[k][i] = c * z[k][i] - s * f;
-						}
-					}
-					if(r == 0.00 && i >= l) continue;
-					d[l] -= p;
-					e[l] = g;
-					e[m] = 0.0;
-				}
-			} while (m != l);
-		}
-		
-		int k;
-		for(int i=0;i<n;i++){
-			p=d[k=i];
-			for(int j=i;j<n;j++){
-				if(d[j] >= p){
-					p=d[k=j];
-				}
-			}
-			if(k!=i){
-				d[k]=d[i];
-				d[i]=p;
-				for(int j=0;j<n;j++){
-					p=z[j][i];
-					z[j][i] = z[j][k];
-					z[j][k] = p;
-				}
-			}
-		}
-	}
-	catch(exception& e) {
-		m->errorOut(e, "PCOACommand", "qtli");
-		exit(1);
-	}
-}
-
-/*********************************************************************************************************************************/
-
 void PCOACommand::output(string fnameRoot, vector<string> name_list, vector<vector<double> >& G, vector<double> d) {
 	try {
 		int rank = name_list.size();