X-Git-Url: https://git.donarmstrong.com/?a=blobdiff_plain;f=pcacommand.cpp;h=2cb37bda8896f8f2cb1f64521f1a73232507f24b;hb=ccae9eef0b44f2d63fdf4a707d0d40243aa1b990;hp=09324c9d3eaf3dd5a525a44cf233abe7df6e7285;hpb=d5bf2c1354d0811a33394d918b15620606560d58;p=mothur.git

diff --git a/pcacommand.cpp b/pcacommand.cpp
index 09324c9..2cb37bd 100644
--- a/pcacommand.cpp
+++ b/pcacommand.cpp
@@ -1,621 +1,468 @@
-
 /*
  *  pcacommand.cpp
- *  Mothur
+ *  mothur
  *
- *  Created by westcott on 1/4/10.
- *  Copyright 2010 Schloss Lab. All rights reserved.
+ *  Created by westcott on 1/7/11.
+ *  Copyright 2011 Schloss Lab. All rights reserved.
  *
  */
 
 #include "pcacommand.h"
+#include "inputdata.h"
 
 //**********************************************************************************************************************
-
-PCACommand::PCACommand(string option){
+vector<string> PCACommand::setParameters(){	
 	try {
-		abort = false;
-		
-		//allow user to run help
-		if(option == "help") { help(); abort = true; }
-		
-		else {
-			//valid paramters for this command
-			string Array[] =  {"phylip"};
-			vector<string> myArray (Array, Array+(sizeof(Array)/sizeof(string)));
-			
-			OptionParser parser(option);
-			map<string, string> parameters = parser. getParameters();
-			
-			ValidParameters validParameter;
-		
-			//check to make sure all parameters are valid for command
-			for (map<string, string>::iterator it = parameters.begin(); it != parameters.end(); it++) { 
-				if (validParameter.isValidParameter(it->first, myArray, it->second) != true) {  abort = true;  }
-			}
-			
-			//required parameters
-			phylipfile = validParameter.validFile(parameters, "phylip", true);
-			if (phylipfile == "not open") { abort = true; }
-			else if (phylipfile == "not found") { phylipfile = ""; abort = true; }	
-			else {	filename = phylipfile;  }
-			
-			//columnfile = validParameter.validFile(parameters, "column", true);
-			//if (columnfile == "not open") { abort = true; }	
-			//else if (columnfile == "not found") { columnfile = ""; }
-			//else {  format = "column";	}
-			
-			//namefile = validParameter.validFile(parameters, "name", true);
-			//if (namefile == "not open") { abort = true; }	
-			//else if (namefile == "not found") { namefile = ""; }
-			
-			
-			//error checking on files	
-			if (phylipfile == "")	{ mothurOut("You must provide a distance file before running the pca command."); mothurOutEndLine(); abort = true; }		
-			//if ((phylipfile == "") && (columnfile == ""))	{ mothurOut("You must provide a distance file before running the pca command."); mothurOutEndLine(); abort = true; }
-			//else if ((phylipfile != "") && (columnfile != "")) { mothurOut("You may not use both the column and the phylip parameters."); mothurOutEndLine(); abort = true; }
-			
-			//if (columnfile != "") {
-			//	if (namefile == "") {  mothurOut("You need to provide a namefile if you are going to use the column format."); mothurOutEndLine(); abort = true; }
-			//}
-			
-			//string temp = validParameter.validFile(parameters, "lt", false);				if (temp == "not found") { temp = "false"; }
-			//bool lt = isTrue(temp);
-			
-			//if (lt)		{  matrix = 2;	}
-			//else		{  matrix = 1;  }
-
-
-		}
-
+		CommandParameter pshared("shared", "InputTypes", "", "", "LRSS", "LRSS", "none",false,false); parameters.push_back(pshared);	
+		CommandParameter prelabund("relabund", "InputTypes", "", "", "LRSS", "LRSS", "none",false,false); parameters.push_back(prelabund);
+		CommandParameter pgroups("groups", "String", "", "", "", "", "",false,false); parameters.push_back(pgroups);
+		CommandParameter pmetric("metric", "Boolean", "", "T", "", "", "",false,false); parameters.push_back(pmetric);
+		CommandParameter plabel("label", "String", "", "", "", "", "",false,false); parameters.push_back(plabel);
+		CommandParameter pinputdir("inputdir", "String", "", "", "", "", "",false,false); parameters.push_back(pinputdir);
+		CommandParameter poutputdir("outputdir", "String", "", "", "", "", "",false,false); parameters.push_back(poutputdir);
+		
+		vector<string> myArray;
+		for (int i = 0; i < parameters.size(); i++) {	myArray.push_back(parameters[i].name);		}
+		return myArray;
 	}
 	catch(exception& e) {
-		errorOut(e, "PCACommand", "PCACommand");
+		m->errorOut(e, "PCACommand", "setParameters");
 		exit(1);
 	}
 }
 //**********************************************************************************************************************
-void PCACommand::help(){
+string PCACommand::getHelpString(){	
 	try {
-	
-		mothurOut("The pca command..."); mothurOutEndLine();
+		string helpString = "";
+		helpString += "The pca command parameters are shared, relabund, label, groups and metric.  shared or relabund is required unless you have a valid current file."; 
+		helpString += "The label parameter is used to analyze specific labels in your input. Default is the first label in your shared or relabund file. Multiple labels may be separated by dashes.\n";
+		helpString += "The groups parameter allows you to specify which groups you would like analyzed. Groupnames are separated by dashes.\n";
+		helpString += "The metric parameter allows you to indicate if would like the pearson correlation coefficient calculated. Default=True";
+		helpString += "Example pca(groups=yourGroups).\n";
+		helpString += "Example pca(groups=A-B-C).\n";
+		helpString += "Note: No spaces between parameter labels (i.e. groups), '=' and parameters (i.e.yourGroups).\n";
+		return helpString;
 	}
 	catch(exception& e) {
-		errorOut(e, "PCACommand", "help");
+		m->errorOut(e, "PCACommand", "getHelpString");
 		exit(1);
 	}
 }
 //**********************************************************************************************************************
-PCACommand::~PCACommand(){}
-//**********************************************************************************************************************
-int PCACommand::execute(){
+string PCACommand::getOutputFileNameTag(string type, string inputName=""){	
 	try {
-	
-		if (abort == true) { return 0; }
-		
-		cout.setf(ios::fixed, ios::floatfield);
-		cout.setf(ios::showpoint);
-		cerr.setf(ios::fixed, ios::floatfield);
-		cerr.setf(ios::showpoint);
-		
-		vector<string> names;
-		vector<vector<double> > D;
-	
-		fbase = filename;
-		if(fbase.find_last_of(".")!=string::npos){
-			fbase.erase(fbase.find_last_of(".")+1); 
-		}
-		else{
-			fbase += ".";
-		}
-		read(filename, names, D);
-   	
-		double offset = 0.0000;
-		vector<double> d;
-		vector<double> e;
-		vector<vector<double> > G = D;
-		vector<vector<double> > copy_G;
-		//int rank = D.size();
-		
-		cout << "\nProcessing...\n";
-		
-		for(int count=0;count<2;count++){
-			recenter(offset, D, G);   
-			tred2(G, d, e);
-			qtli(d, e, G);
-			offset = d[d.size()-1];
-			if(offset > 0.0) break;
-		} 
-		
-		
-		output(fbase, names, G, d);
-		
-		return 0;
+        string outputFileName = "";
+		map<string, vector<string> >::iterator it;
+        
+        //is this a type this command creates
+        it = outputTypes.find(type);
+        if (it == outputTypes.end()) {  m->mothurOut("[ERROR]: this command doesn't create a " + type + " output file.\n"); }
+        else {
+            if (type == "pca") {  outputFileName =  "pca.axes"; }
+            else if (type == "loadings") {  outputFileName =  "pca.loadings"; }
+            else { m->mothurOut("[ERROR]: No definition for type " + type + " output file tag.\n"); m->control_pressed = true;  }
+        }
+        return outputFileName;
 	}
 	catch(exception& e) {
-		errorOut(e, "PCACommand", "execute");
+		m->errorOut(e, "PCACommand", "getOutputFileNameTag");
 		exit(1);
 	}
 }
-/*********************************************************************************************************************************/
 
-inline double SIGN(const double a, const double b)
-{
-    return b>=0 ? (a>=0 ? a:-a) : (a>=0 ? -a:a);
-}
-
-/*********************************************************************************************************************************/
-
-void PCACommand::get_comment(istream& f, char begin, char end){
+//**********************************************************************************************************************
+PCACommand::PCACommand(){	
 	try {
-		char d=f.get();
-		while(d != end){	d = f.get();	}
-		d = f.peek();
+		abort = true; calledHelp = true; 
+		setParameters();
+		vector<string> tempOutNames;
+		outputTypes["pca"] = tempOutNames;
+		outputTypes["loadings"] = tempOutNames;
 	}
 	catch(exception& e) {
-		errorOut(e, "PCACommand", "get_comment");
+		m->errorOut(e, "PCACommand", "PCACommand");
 		exit(1);
 	}
-}	
-
-/*********************************************************************************************************************************/
+}
+//**********************************************************************************************************************
 
-void PCACommand::read_mega(istream& f, vector<string>& name_list, vector<vector<double> >& d){
+PCACommand::PCACommand(string option)  {
 	try {
-		get_comment(f, '#', '\n');
+		abort = false; calledHelp = false;   
 		
-		char test = f.peek();
+		//allow user to run help
+		if(option == "help") { help(); abort = true; calledHelp = true; }
+		else if(option == "citation") { citation(); abort = true; calledHelp = true;}
 		
-		while(test == '!'){								//get header comments
-			get_comment(f, '!', ';');
-			while(isspace(test=f.get()))		{;}
-			f.putback(test);
-			test = f.peek();
-		}
-		while(test != '\n'){							//get sequence names
-			get_comment(f, '[', ']');
-			char d = f.get();
-			d = f.get();
-			if(d == '#'){
-				string name;
-				f >> name;
-				name_list.push_back(name);
-				while(isspace(test=f.get()))		{;}
-				f.putback(test);
-			}
-			else{
-				break;
-			}
-		}
-		int rank = name_list.size();
-		d.resize(rank);
-		for(int i=0;i<rank;i++){		d[i].resize(rank);		}
-		
-		d[0][0] = 0.0000;
-		get_comment(f, '[', ']');
-		for(int i=1;i<rank;i++){
-			get_comment(f, '[', ']');
-			d[i][i]=0.0000;
-			for(int j=0;j<i;j++){
-				f >> d[i][j];
-				if (d[i][j] == -0.0000)
-					d[i][j] = 0.0000;
-				d[j][i]=d[i][j];
+		else {
+			vector<string> myArray = setParameters();
+			
+			OptionParser parser(option);
+			map<string, string> parameters = parser. getParameters();
+			
+			ValidParameters validParameter;
+			map<string, string>::iterator it;
+			
+			//check to make sure all parameters are valid for command
+			for (it = parameters.begin(); it != parameters.end(); it++) { 
+				if (validParameter.isValidParameter(it->first, myArray, it->second) != true) {  abort = true;  }
 			}
-		}
-	}
-	catch(exception& e) {
-		errorOut(e, "PCACommand", "read_mega");
-		exit(1);
-	}
-}
-
-/*********************************************************************************************************************************/
-
-void PCACommand::read_phylip(istream& f, int square_m, vector<string>& name_list, vector<vector<double> >& d){
-	try {
-		//     int count1=0;
-		//     int count2=0;
-		
-		int rank;
-		f >> rank;
-		
-		name_list.resize(rank);
-		d.resize(rank);
-		if(square_m == 1){
-			for(int i=0;i<rank;i++)
-				d[i].resize(rank);
-			for(int i=0;i<rank;i++) {
-				f >> name_list[i];
-				//			cout << i << "\t" << name_list[i] << endl;
-				for(int j=0;j<rank;j++) {
-					f >> d[i][j];
-					if (d[i][j] == -0.0000)
-						d[i][j] = 0.0000;
+	
+			//initialize outputTypes
+			vector<string> tempOutNames;
+			outputTypes["pca"] = tempOutNames;
+			outputTypes["loadings"] = tempOutNames;
+			
+			//if the user changes the input directory command factory will send this info to us in the output parameter 
+			string inputDir = validParameter.validFile(parameters, "inputdir", false);		
+			if (inputDir == "not found"){	inputDir = "";		}
+			else {
+				string path;
+				it = parameters.find("shared");
+				//user has given a template file
+				if(it != parameters.end()){ 
+					path = m->hasPath(it->second);
+					//if the user has not given a path then, add inputdir. else leave path alone.
+					if (path == "") {	parameters["shared"] = inputDir + it->second;		}
+				}
+				
+				it = parameters.find("relabund");
+				//user has given a template file
+				if(it != parameters.end()){ 
+					path = m->hasPath(it->second);
+					//if the user has not given a path then, add inputdir. else leave path alone.
+					if (path == "") {	parameters["relabund"] = inputDir + it->second;		}
 				}
 			}
-		}
-		else if(square_m == 2){
-			for(int i=0;i<rank;i++){
-				d[i].resize(rank);
-			}
-			d[0][0] = 0.0000;
-			f >> name_list[0];
-			for(int i=1;i<rank;i++){
-				f >> name_list[i];
-				d[i][i]=0.0000;
-				for(int j=0;j<i;j++){
-					f >> d[i][j];
-					if (d[i][j] == -0.0000)
-						d[i][j] = 0.0000;
-					d[j][i]=d[i][j];
+			
+			//check for required parameters
+			sharedfile = validParameter.validFile(parameters, "shared", true);
+			if (sharedfile == "not open") { sharedfile = ""; abort = true; }	
+			else if (sharedfile == "not found") { sharedfile = ""; }
+			else {  mode = "sharedfile"; inputFile = sharedfile; m->setSharedFile(sharedfile); }
+			
+			relabundfile = validParameter.validFile(parameters, "relabund", true);
+			if (relabundfile == "not open") { relabundfile = ""; abort = true; }	
+			else if (relabundfile == "not found") { relabundfile = ""; }
+			else {  mode = "relabund"; inputFile = relabundfile; m->setRelAbundFile(relabundfile); }
+			
+			
+			if ((sharedfile == "") && (relabundfile == "")) { 
+				//is there are current file available for any of these?
+				//give priority to shared, then list, then rabund, then sabund
+				//if there is a current shared file, use it
+				sharedfile = m->getSharedFile(); 
+				if (sharedfile != "") { inputFile = sharedfile; mode = "sharedfile"; m->mothurOut("Using " + sharedfile + " as input file for the shared parameter."); m->mothurOutEndLine(); }
+				else { 
+					relabundfile = m->getRelAbundFile(); 
+					if (relabundfile != "") { inputFile = relabundfile; mode = "relabund"; m->mothurOut("Using " + relabundfile + " as input file for the relabund parameter."); m->mothurOutEndLine(); }
+					else { 
+						m->mothurOut("No valid current files. You must provide a relabund or shared file."); m->mothurOutEndLine(); 
+						abort = true;
+					}
 				}
 			}
+				
+			//if the user changes the output directory command factory will send this info to us in the output parameter 
+			outputDir = validParameter.validFile(parameters, "outputdir", false);		if (outputDir == "not found"){	
+				outputDir = "";	
+				outputDir += m->hasPath(inputFile); //if user entered a file with a path then preserve it	
+			}
+						
+			string temp = validParameter.validFile(parameters, "metric", false);	if (temp == "not found"){	temp = "T";				}
+			metric = m->isTrue(temp); 
+			
+			label = validParameter.validFile(parameters, "label", false);			
+			if (label == "not found") { label = ""; if(labels.size() == 0) {  m->mothurOut("You did not provide a label, I will use the first label in your inputfile."); m->mothurOutEndLine(); } }
+			else { m->splitAtDash(label, labels); }
+			
+			groups = validParameter.validFile(parameters, "groups", false);			
+			if (groups == "not found") { groups = "";  }
+			else { m->splitAtDash(groups, Groups);	}			
+			m->setGroups(Groups);			
+			
 		}
+		
 	}
 	catch(exception& e) {
-		errorOut(e, "PCACommand", "read_phylip");
+		m->errorOut(e, "PCACommand", "PCACommand");
 		exit(1);
 	}
-
 }
-
-/*********************************************************************************************************************************/
-
-void PCACommand::read(string fname, vector<string>& names, vector<vector<double> >& D){
+//**********************************************************************************************************************
+int PCACommand::execute(){
 	try {
-		ifstream f;
-		openInputFile(fname, f);
 		
-		char test = f.peek();
+		if (abort == true) { if (calledHelp) { return 0; }  return 2;	}
+		
+		cout.setf(ios::fixed, ios::floatfield);
+		cout.setf(ios::showpoint);
+		cerr.setf(ios::fixed, ios::floatfield);
+		cerr.setf(ios::showpoint);
+		
+		//get first line of shared file
+		vector< vector<double> > matrix;
+		InputData* input;
+		if (mode == "sharedfile")			{  
+			input = new InputData(inputFile, "sharedfile");
+		}else if (mode == "relabund")	{ 
+			input = new InputData(inputFile, "relabund");
+		}else {  m->mothurOut("[ERROR]: filetype not recognized."); m->mothurOutEndLine();  return 0; }
+		
+		vector<SharedRAbundFloatVector*> lookupFloat = input->getSharedRAbundFloatVectors();
+		string lastLabel = lookupFloat[0]->getLabel();
+			
+		set<string> processedLabels;
+		set<string> userLabels = labels;
 		
-		if(test == '#'){
-			read_mega(f, names, D);
+		//if the user gave no labels, then use the first one read
+		if (labels.size() == 0) { 
+			label = lastLabel;  
+			
+			process(lookupFloat);
 		}
-		else{
-			//check whether matrix is square
-			char d;
-			int m = 1;
-			int numSeqs;
-			string name;
+		
+		//as long as you are not at the end of the file or done wih the lines you want
+		while((lookupFloat[0] != NULL) && (userLabels.size() != 0)) {
 			
-			f >> numSeqs >> name; 
+			if (m->control_pressed) {  for (int i = 0; i < outputNames.size(); i++) {	m->mothurRemove(outputNames[i]);  } delete input; for (int i = 0; i < lookupFloat.size(); i++) {  delete lookupFloat[i];  }  lookupFloat.clear(); return 0;  }
 			
-			while((d=f.get()) != EOF){
+			if(labels.count(lookupFloat[0]->getLabel()) == 1){
+				processedLabels.insert(lookupFloat[0]->getLabel());
+				userLabels.erase(lookupFloat[0]->getLabel());
 				
-				//is d a number meaning its square
-				if(isalnum(d)){ 
-					m = 1; 
-					break; 
-				}
+				process(lookupFloat);
+			}
+			
+			if ((m->anyLabelsToProcess(lookupFloat[0]->getLabel(), userLabels, "") == true) && (processedLabels.count(lastLabel) != 1)) {
+				string saveLabel = lookupFloat[0]->getLabel();
 				
-				//is d a line return meaning its lower triangle
-				if(d == '\n'){
-					m = 2;
-					break;
-				}
+				for (int i = 0; i < lookupFloat.size(); i++) {  delete lookupFloat[i];  }  lookupFloat.clear();
+				lookupFloat = input->getSharedRAbundFloatVectors(lastLabel);
+				
+				process(lookupFloat);
+				
+				processedLabels.insert(lookupFloat[0]->getLabel());
+				userLabels.erase(lookupFloat[0]->getLabel());
+				
+				//restore real lastlabel to save below
+				lookupFloat[0]->setLabel(saveLabel);
 			}
-			f.close();
 			
-			//reopen to get back to beginning
-			openInputFile(fname, f);			
-			read_phylip(f, m, names, D);
+			lastLabel = lookupFloat[0]->getLabel();			
+			
+			//get next line to process
+			//prevent memory leak
+			for (int i = 0; i < lookupFloat.size(); i++) {  delete lookupFloat[i];  } lookupFloat.clear();
+			lookupFloat = input->getSharedRAbundFloatVectors();
 		}
 		
-		//int rank = D.size();
-	}
-	catch(exception& e) {
-		errorOut(e, "PCACommand", "read");
-		exit(1);
-	}
-}
-
-/*********************************************************************************************************************************/
-double PCACommand::pythag(double a, double b){
-    return(pow(a*a+b*b,0.5));
-}
-/*********************************************************************************************************************************/
-
-void PCACommand::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);
-		}
+		if (m->control_pressed) { for (int i = 0; i < outputNames.size(); i++) {	m->mothurRemove(outputNames[i]);  } delete input; for (int i = 0; i < lookupFloat.size(); i++) {  delete lookupFloat[i];  } lookupFloat.clear(); return 0;  }
 		
-		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];
-				}
+		//output error messages about any remaining user labels
+		set<string>::iterator it;
+		bool needToRun = false;
+		for (it = userLabels.begin(); it != userLabels.end(); it++) {  
+			m->mothurOut("Your file does not include the label " + *it); 
+			if (processedLabels.count(lastLabel) != 1) {
+				m->mothurOut(". I will use " + lastLabel + "."); m->mothurOutEndLine();
+				needToRun = true;
+			}else {
+				m->mothurOut(". Please refer to " + lastLabel + "."); m->mothurOutEndLine();
 			}
 		}
+		
+		//run last label if you need to
+		if (needToRun == true)  {
+			for (int i = 0; i < lookupFloat.size(); i++) {  if (lookupFloat[i] != NULL) {	delete lookupFloat[i];	} }  lookupFloat.clear();
+			lookupFloat = input->getSharedRAbundFloatVectors(lastLabel);
+			
+			process(lookupFloat);
+			
+			for (int i = 0; i < lookupFloat.size(); i++) {  if (lookupFloat[i] != NULL) {	delete lookupFloat[i];	} } lookupFloat.clear();
+		}	
+		
+		for (int i = 0; i < lookupFloat.size(); i++) {  if (lookupFloat[i] != NULL) {	delete lookupFloat[i];	} } lookupFloat.clear();
+		delete input;
+		
+		if (m->control_pressed) { for (int i = 0; i < outputNames.size(); i++) {	m->mothurRemove(outputNames[i]);  } return 0; }
+		
+		m->mothurOutEndLine();
+		m->mothurOut("Output File Names: "); m->mothurOutEndLine();
+		for (int i = 0; i < outputNames.size(); i++) {	m->mothurOut(outputNames[i]); m->mothurOutEndLine();	}
+		m->mothurOutEndLine();
+		
+		return 0;
 	}
 	catch(exception& e) {
-		errorOut(e, "PCACommand", "matrix_mult");
+		m->errorOut(e, "PCACommand", "execute");
 		exit(1);
 	}
-
 }
 
-/*********************************************************************************************************************************/
-
-void PCACommand::recenter(double offset, vector<vector<double> > D, vector<vector<double> >& G){
+/**********************************************************************************************************************
+vector< vector<double> > PCACommand::createMatrix(vector<SharedRAbundFloatVector*> lookupFloat){
 	try {
-		int rank = D.size();
+		vector< vector<double> > matrix; matrix.resize(lookupFloat.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);
+		//fill matrix with shared files relative abundances
+		for (int i = 0; i < lookupFloat.size(); i++) {
+			for (int j = 0; j < lookupFloat[i]->getNumBins(); j++) {
+				matrix[i].push_back(lookupFloat[i]->getAbundance(j));
+			}
 		}
 		
-		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;
+		vector< vector<double> > transposeMatrix; transposeMatrix.resize(matrix[0].size());
+		for (int i = 0; i < transposeMatrix.size(); i++) {
+			for (int j = 0; j < matrix.size(); j++) {
+				transposeMatrix[i].push_back(matrix[j][i]);
 			}
 		}
 		
-		matrix_mult(C,A,A);
-		matrix_mult(A,C,G);
+		matrix = linearCalc.matrix_mult(matrix, transposeMatrix);
+		
+		return matrix;
 	}
 	catch(exception& e) {
-		errorOut(e, "PCACommand", "recenter");
+		m->errorOut(e, "PCACommand", "createMatrix");	
 		exit(1);
 	}
+}*/
+//**********************************************************************************************************************
 
-}
-
-/*********************************************************************************************************************************/
-
-//  This function is taken from Numerical Recipes in C++ by Press et al., 2nd edition, pg. 479
-
-void PCACommand::tred2(vector<vector<double> >& a, vector<double>& d, vector<double>& e){
+int PCACommand::process(vector<SharedRAbundFloatVector*>& lookupFloat){
 	try {
-		double scale, hh, h, g, f;
-		
-		int n = a.size();
+		m->mothurOut("\nProcessing " + lookupFloat[0]->getLabel()); m->mothurOutEndLine();
+	
+		int numOTUs = lookupFloat[0]->getNumBins();
+		int numSamples = lookupFloat.size();
 		
-		d.resize(n);
-		e.resize(n);
+		vector< vector<double> > matrix(numSamples);
+		vector<double> colMeans(numOTUs);
 		
-		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];
-			}
+		//fill matrix with shared relative abundances, re-center
+		for (int i = 0; i < lookupFloat.size(); i++) {
+			matrix[i].resize(numOTUs, 0);
 			
-			d[i] = h;
+			for (int j = 0; j < numOTUs; j++) {
+				matrix[i][j] = lookupFloat[i]->getAbundance(j);
+				colMeans[j] += matrix[i][j];
+			}
 		}
 		
-		d[0] = 0.0000;
-		e[0] = 0.0000;
+
+		for(int j=0;j<numOTUs;j++){
+			colMeans[j] = colMeans[j] / (double)numSamples;
+		}
 		
-		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;
+		vector<vector<double> > centered = matrix;
+		for(int i=0;i<numSamples;i++){
+			for(int j=0;j<numOTUs;j++){
+				centered[i][j] = centered[i][j] - colMeans[j];				
 			}
 		}
-	}
-	catch(exception& e) {
-		errorOut(e, "PCACommand", "tred2");
-		exit(1);
-	}
-
-}
 
-/*********************************************************************************************************************************/
+		
+		vector< vector<double> > transpose(numOTUs);
+		for (int i = 0; i < numOTUs; i++) {
+			transpose[i].resize(numSamples, 0);
+			
+			for (int j = 0; j < numSamples; j++) {
+				transpose[i][j] = centered[j][i];
+			}
+		}
 
-//  This function is taken from Numerical Recipes in C++ by Press et al., 2nd edition, pg. 479
+		vector<vector<double> > crossProduct = linearCalc.matrix_mult(transpose, centered);	
+		
+		vector<double> d;
+		vector<double> e;
 
-void PCACommand::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;
+		linearCalc.tred2(crossProduct, d, e);		if (m->control_pressed) { return 0; }
+		linearCalc.qtli(d, e, crossProduct);		if (m->control_pressed) { return 0; }
 		
-		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);
-		}
+		vector<vector<double> > X = linearCalc.matrix_mult(centered, crossProduct);
 		
-		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;
-				}
+		if (m->control_pressed) { return 0; }
+		
+		string fbase = outputDir + m->getRootName(m->getSimpleName(inputFile));
+		string outputFileName = fbase + lookupFloat[0]->getLabel();
+		output(outputFileName, m->getGroups(), X, d);
+		
+		if (metric) {   
+			
+			vector<vector<double> > observedEuclideanDistance = linearCalc.getObservedEuclideanDistance(centered);
+			
+			for (int i = 1; i < 4; i++) {
+				
+				vector< vector<double> > PCAEuclidDists = linearCalc.calculateEuclidianDistance(X, i); //G is the pca file
+				
+				if (m->control_pressed) { for (int i = 0; i < outputNames.size(); i++) {	m->mothurRemove(outputNames[i]);  } return 0; }
+
+				double corr = linearCalc.calcPearson(PCAEuclidDists, observedEuclideanDistance);
+								
+				m->mothurOut("Rsq " + toString(i) + " axis: " + toString(corr * corr)); m->mothurOutEndLine();
+				
+				if (m->control_pressed) { for (int i = 0; i < outputNames.size(); i++) {	m->mothurRemove(outputNames[i]);  } return 0; }
 			}
 		}
+		
+		return 0;
 	}
 	catch(exception& e) {
-		errorOut(e, "PCACommand", "qtli");
+		m->errorOut(e, "PCACommand", "process");	
 		exit(1);
 	}
 }
-
 /*********************************************************************************************************************************/
 
-void PCACommand::output(string fnameRoot, vector<string> name_list, vector<vector<double> > G, vector<double> d) {
+void PCACommand::output(string fnameRoot, vector<string> name_list, vector<vector<double> >& G, vector<double> d) {
 	try {
-		int rank = name_list.size();
+
+		int numEigenValues = d.size();
 		double dsum = 0.0000;
-		for(int i=0;i<rank;i++){
+		for(int i=0;i<numEigenValues;i++){
 			dsum += d[i];
-			for(int j=0;j<rank;j++){
-				if(d[j] >= 0)	{	G[i][j] *= pow(d[j],0.5);	}
-				else			{	G[i][j] = 0.00000;			}
-			}
 		}
 		
-		ofstream pcaData((fnameRoot+"pca").c_str(), ios::trunc);
+		ofstream pcaData;
+        m->openOutputFile((fnameRoot+"."+getOutputFileNameTag("pca")), pcaData);
 		pcaData.setf(ios::fixed, ios::floatfield);
 		pcaData.setf(ios::showpoint);	
+		outputNames.push_back(fnameRoot+"."+getOutputFileNameTag("pca"));
+		outputTypes["pca"].push_back(fnameRoot+"."+getOutputFileNameTag("pca"));
 		
-		ofstream pcaLoadings((fnameRoot+"pca.loadings").c_str(), ios::trunc);
+		ofstream pcaLoadings;
+         m->openOutputFile((fnameRoot+"."+getOutputFileNameTag("loadings")), pcaLoadings);
 		pcaLoadings.setf(ios::fixed, ios::floatfield);
-		pcaLoadings.setf(ios::showpoint);	
+		pcaLoadings.setf(ios::showpoint);
+		outputNames.push_back(fnameRoot+"."+getOutputFileNameTag("loadings"));
+		outputTypes["loadings"].push_back(fnameRoot+"."+getOutputFileNameTag("loadings"));	
 		
 		pcaLoadings << "axis\tloading\n";
-		for(int i=0;i<rank;i++){
+		for(int i=0;i<numEigenValues;i++){
 			pcaLoadings << i+1 << '\t' << d[i] * 100.0 / dsum << endl;
 		}
 		
-		pcaData << "SeqName";
-		for(int i=0;i<rank;i++){
+		pcaData << "group";
+		for(int i=0;i<numEigenValues;i++){
 			pcaData << '\t' << "axis" << i+1;
 		}
 		pcaData << endl;
 		
-		for(int i=0;i<rank;i++){
+		for(int i=0;i<name_list.size();i++){
 			pcaData << name_list[i] << '\t';
-			for(int j=0;j<rank;j++){
+			for(int j=0;j<numEigenValues;j++){
 				pcaData << G[i][j] << '\t';
 			}
 			pcaData << endl;
 		}
 	}
 	catch(exception& e) {
-		errorOut(e, "PCACommand", "output");
+		m->errorOut(e, "PCACommand", "output");
 		exit(1);
 	}
 }
-
 /*********************************************************************************************************************************/
 
-void PCACommand::print_matrix(vector<vector<double> > A) {
-	try {
-		int rank = A.size();
-		for(int i=0;i<rank;i++){
-			for(int j=0;j<rank;j++){
-				cout << A[i][j] << "  ";
-			}
-			cout << endl;
-		}
-	}
-	catch(exception& e) {
-		errorOut(e, "PCACommand", "print_matrix");
-		exit(1);
-	}
-}
-/*********************************************************************************************************************************/