X-Git-Url: https://git.donarmstrong.com/?a=blobdiff_plain;f=distancecommand.cpp;h=065c14a1671a8abe78e1dda115208936185f496e;hb=cd37904452dc95b183ff313ff05720c562902487;hp=6fee651cb8963680ec294051e671cfc57b9c6135;hpb=cdcf99b7760701e6869ca8e4e6e91c8e8c4ae186;p=mothur.git

diff --git a/distancecommand.cpp b/distancecommand.cpp
index 6fee651..065c14a 100644
--- a/distancecommand.cpp
+++ b/distancecommand.cpp
@@ -16,42 +16,83 @@
 
 //**********************************************************************************************************************
 
-DistanceCommand::DistanceCommand(){
+DistanceCommand::DistanceCommand(string option){
 	try {
-		globaldata = GlobalData::getInstance();
-		validCalculator = new ValidCalculators();
-		countends = globaldata->getCountEnds();
-		convert(globaldata->getProcessors(), processors);
-		convert(globaldata->getCutOff(), cutoff);
+		abort = false;
+		Estimators.clear();
 		
-		int i;
-		if (countends == "T") {
-			for (i=0; i<globaldata->Estimators.size(); i++) {
-				if (validCalculator->isValidCalculator("distance", globaldata->Estimators[i]) == true) { 
-					if (globaldata->Estimators[i] == "nogaps") { 
-						distCalculator = new ignoreGaps();
-					}else if (globaldata->Estimators[i] == "eachgap") { 
-						distCalculator = new eachGapDist();	
-					}else if (globaldata->Estimators[i] == "onegap") {
-						distCalculator = new oneGapDist();					}
-				}
+		//allow user to run help
+		if(option == "help") { help(); abort = true; }
+		
+		else {
+			//valid paramters for this command
+			string Array[] =  {"fasta", "phylip", "calc", "countends", "cutoff", "processors"};
+			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 it2 = parameters.begin(); it2 != parameters.end(); it2++) { 
+				if (validParameter.isValidParameter(it2->first, myArray, it2->second) != true) {  abort = true;  }
+			}
+			
+			//check for required parameters
+			fastafile = validParameter.validFile(parameters, "fasta", true);
+			if (fastafile == "not found") { cout << "fasta is a required parameter for the dist.seqs command." << endl; abort = true; }
+			else if (fastafile == "not open") { abort = true; }	
+			else{
+				ifstream inFASTA;
+				openInputFile(fastafile, inFASTA);
+				alignDB = SequenceDB(inFASTA); 
 			}
-		}else {
-			for (i=0; i<globaldata->Estimators.size(); i++) {
-				if (validCalculator->isValidCalculator("distance", globaldata->Estimators[i]) == true) { 
-					if (globaldata->Estimators[i] == "nogaps") { 
-						distCalculator = new ignoreGaps();	
-					}else if (globaldata->Estimators[i] == "eachgap") { 
-						distCalculator = new eachGapIgnoreTermGapDist();
-					}else if (globaldata->Estimators[i] == "onegap") { 
-						distCalculator = new oneGapIgnoreTermGapDist();	
+
+			//check for optional parameter and set defaults
+			// ...at some point should added some additional type checking...
+			calc = validParameter.validFile(parameters, "calc", false);			
+			if (calc == "not found") { calc = "onegap";  }
+			else { 
+				 if (calc == "default")  {  calc = "onegap";  }
+			}
+			splitAtDash(calc, Estimators);
+
+			string temp;
+			temp = validParameter.validFile(parameters, "countends", false);	if(temp == "not found"){	temp = "T";	}
+			convert(temp, countends); 
+			
+			temp = validParameter.validFile(parameters, "cutoff", false);		if(temp == "not found"){	temp = "1.0"; }
+			convert(temp, cutoff); 
+			
+			temp = validParameter.validFile(parameters, "processors", false);	if(temp == "not found"){	temp = "1"; }
+			convert(temp, processors); 
+			
+			phylip = validParameter.validFile(parameters, "phylip", false);		if(phylip == "not found"){	phylip = "F"; }
+	
+			
+			ValidCalculators validCalculator;
+			
+			if (isTrue(countends) == true) {
+				for (int i=0; i<Estimators.size(); i++) {
+					if (validCalculator.isValidCalculator("distance", Estimators[i]) == true) { 
+						if (Estimators[i] == "nogaps")			{	distCalculator = new ignoreGaps();	}
+						else if (Estimators[i] == "eachgap")	{	distCalculator = new eachGapDist();	}
+						else if (Estimators[i] == "onegap")		{	distCalculator = new oneGapDist();	}
+					}
+				}
+			}else {
+				for (int i=0; i<Estimators.size(); i++) {
+					if (validCalculator.isValidCalculator("distance", Estimators[i]) == true) { 
+						if (Estimators[i] == "nogaps")		{	distCalculator = new ignoreGaps();					}
+						else if (Estimators[i] == "eachgap"){	distCalculator = new eachGapIgnoreTermGapDist();	}
+						else if (Estimators[i] == "onegap")	{	distCalculator = new oneGapIgnoreTermGapDist();		}
 					}
 				}
 			}
+
 		}
 				
-		//reset calc for next command
-		globaldata->setCalc("");
 	}
 	catch(exception& e) {
 		cout << "Standard Error: " << e.what() << " has occurred in the DistanceCommand class Function DistanceCommand. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
@@ -62,122 +103,106 @@ DistanceCommand::DistanceCommand(){
 		exit(1);
 	}	
 }
+
+//**********************************************************************************************************************
+
+DistanceCommand::~DistanceCommand(){
+	
+	for(int i=0;i<lines.size();i++){
+		delete lines[i];
+	}
+	
+}
+	
+//**********************************************************************************************************************
+
+void DistanceCommand::help(){
+	try {
+		cout << "The dist.seqs command reads a file containing sequences and creates a distance file." << "\n";
+		cout << "The dist.seqs command parameters are fasta, calc, countends, cutoff and processors.  " << "\n";
+		cout << "The fasta parameter is required." << "\n";
+		cout << "The calc parameter allows you to specify the method of calculating the distances.  Your options are: nogaps, onegap or eachgap. The default is onegap." << "\n";
+		cout << "The countends parameter allows you to specify whether to include terminal gaps in distance.  Your options are: T or F. The default is T." << "\n";
+		cout << "The cutoff parameter allows you to specify maximum distance to keep. The default is 1.0." << "\n";
+		cout << "The processors parameter allows you to specify number of processors to use.  The default is 1." << "\n";
+		cout << "The dist.seqs command should be in the following format: " << "\n";
+		cout << "dist.seqs(fasta=yourFastaFile, calc=yourCalc, countends=yourEnds, cutoff= yourCutOff, processors=yourProcessors) " << "\n";
+		cout << "Example dist.seqs(fasta=amazon.fasta, calc=eachgap, countends=F, cutoff= 2.0, processors=3)." << "\n";
+		cout << "Note: No spaces between parameter labels (i.e. calc), '=' and parameters (i.e.yourCalc)." << "\n" << "\n";
+	}
+	catch(exception& e) {
+		cout << "Standard Error: " << e.what() << " has occurred in the DistanceCommand class Function help. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
+		exit(1);
+	}
+	catch(...) {
+		cout << "An unknown error has occurred in the DistanceCommand class function help. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
+		exit(1);
+	}	
+}
+
+
 //**********************************************************************************************************************
 
 int DistanceCommand::execute(){
 	try {
 		
-		//read file
-		string filename = globaldata->inputFileName;
+		if (abort == true) { return 0; }
 		
-		if(globaldata->getFastaFile() != "") {
-			readSeqs =  new ReadFasta(filename); }
-		else if(globaldata->getNexusFile() != "") {
-			readSeqs = new ReadNexus(filename); }
-		else if(globaldata->getClustalFile() != "") {
-			readSeqs = new ReadClustal(filename); }
-		else if(globaldata->getPhylipFile() != "") {
-			readSeqs = new ReadPhylip(filename); }
-			
-		readSeqs->read();
-		seqDB = readSeqs->getDB();
-	
-		int numSeqs = seqDB->getNumSeqs();
+		int numSeqs = alignDB.getNumSeqs();
 		cutoff += 0.005;
-
-		string distFile = getRootName(globaldata->getFastaFile()) + "dist";
 		
-		remove(distFile.c_str());
+		string outputFile;
 		
+		//doses the user want the phylip formatted file as well
+		if (isTrue(phylip) == true) {
+			outputFile = getRootName(fastafile) + "phylip.dist";
+			remove(outputFile.c_str());
+			
+			//output numSeqs to phylip formatted dist file
+		}else { //user wants column format
+			outputFile = getRootName(fastafile) + "dist";
+			remove(outputFile.c_str());
+		}
+				
 		//#	if defined (_WIN32)
-			//figure out how to implement the fork and wait commands in windows
-		//	driver(distCalculator, seqDB, 0, numSeqs, distFile, cutoff);
+		//figure out how to implement the fork and wait commands in windows
+		//	driver(distCalculator, seqDB, 0, numSeqs, distFile, phylipFile, cutoff);
 		//#	endif
 		
-		#	if defined (__APPLE__) || (__MACH__) || (linux) || (__linux)
-			if(processors == 1){
-				driver(distCalculator, seqDB, 0, numSeqs, distFile, cutoff);	
-			}	
-			else if(processors == 2){
-		
-				int pid = fork();
-				if(pid > 0){
-					driver(distCalculator, seqDB, 0, (numSeqs/sqrt(2)), distFile + "tempa", cutoff);	
-					appendFiles((distFile+"tempa"), distFile);
-					remove((distFile + "tempa").c_str());	
-				}
-				else{
-					driver(distCalculator, seqDB, (numSeqs/sqrt(2)), numSeqs, distFile + "tempb", cutoff);	
-					appendFiles((distFile+"tempb"), distFile);
-					remove((distFile + "tempb").c_str());	
-				}
-				wait(NULL);
-
-			}
-			else if(processors == 3){
-				int pid1 = fork();
-				if(pid1 > 0){
-					int pid2 = fork();
-					if(pid2 > 0){
-						driver(distCalculator, seqDB, 0, sqrt(3) * numSeqs / 3, distFile + "tempa", cutoff);
-						appendFiles(distFile+"tempa", distFile);
-						remove((distFile + "tempa").c_str());	
-					}
-					else{
-						driver(distCalculator, seqDB, sqrt(3) * numSeqs / 3, sqrt(6) * numSeqs / 3, distFile + "tempb", cutoff);	
-						appendFiles(distFile+"tempb", distFile);
-						remove((distFile + "tempb").c_str());				
-					}
-					wait(NULL);
-				}
-				else{
-					driver(distCalculator, seqDB, sqrt(6) * numSeqs / 3, numSeqs, distFile + "tempc", cutoff);	
-					appendFiles(distFile+"tempc", distFile);
-					remove((distFile + "tempc").c_str());			
-				}
-				wait(NULL);
+				
+#if defined (__APPLE__) || (__MACH__) || (linux) || (__linux)
+		//if you don't need to fork anything
+		if(processors == 1){
+			driver(0, numSeqs, outputFile, cutoff);
+		}else{ //you have multiple processors
+			
+			for (int i = 0; i < processors; i++) {
+				lines.push_back(new linePair());
+				lines[i]->start = int (sqrt(float(i)/float(processors)) * numSeqs);
+				lines[i]->end = int (sqrt(float(i+1)/float(processors)) * numSeqs);
 			}
-			else if(processors == 4){
-				int pid1 = fork();
-				if(pid1 > 0){
-					int pid2 = fork();
-					if(pid2 > 0){
-						driver(distCalculator, seqDB, 0, numSeqs / 2, distFile + "tempa", cutoff);	
-						appendFiles(distFile+"tempa", distFile);
-						remove((distFile + "tempa").c_str());			
-					}
-					else{
-						driver(distCalculator, seqDB, numSeqs / 2, (numSeqs/sqrt(2)), distFile + "tempb", cutoff);	
-						appendFiles(distFile+"tempb", distFile);
-						remove((distFile + "tempb").c_str());				
-					}
-					wait(NULL);
-				}
-				else{
-					int pid3 = fork();
-					if(pid3 > 0){
-						driver(distCalculator, seqDB, (numSeqs/sqrt(2)), (sqrt(3) * numSeqs / 2), distFile + "tempc", cutoff);	
-						appendFiles(distFile+"tempc", distFile);
-						remove((distFile + "tempc").c_str());				
-					}
-					else{
-						driver(distCalculator, seqDB, (sqrt(3) * numSeqs / 2), numSeqs, distFile + "tempd", cutoff);	
-						appendFiles(distFile+"tempd", distFile);
-						remove((distFile + "tempd").c_str());				
-					}
-					wait(NULL);
-				}
-				wait(NULL);
+
+			createProcesses(outputFile); 
+		
+			map<int, int>::iterator it = processIDS.begin();
+			rename((outputFile + toString(it->second) + ".temp").c_str(), outputFile.c_str());
+			it++;
+			
+			//append and remove temp files
+			for (; it != processIDS.end(); it++) {
+				appendFiles((outputFile + toString(it->second) + ".temp"), outputFile);
+				remove((outputFile + toString(it->second) + ".temp").c_str());
 			}
-			wait(NULL);
-		#	else
-			driver(distCalculator, seqDB, 0, numSeqs, distFile, cutoff);
-		#	endif
-	
+		}
+#else
+		ifstream inFASTA
+		driver(0, numSeqs, outputFile, cutoff);
+#endif
+		
 		delete distCalculator;
-	
+		
 		return 0;
-
+		
 	}
 	catch(exception& e) {
 		cout << "Standard Error: " << e.what() << " has occurred in the DistanceCommand class Function execute. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
@@ -188,35 +213,80 @@ int DistanceCommand::execute(){
 		exit(1);
 	}	
 }
+/**************************************************************************************************/
+void DistanceCommand::createProcesses(string filename) {
+	try {
+#if defined (__APPLE__) || (__MACH__) || (linux) || (__linux)
+		int process = 0;
+		processIDS.clear();
+		
+		//loop through and create all the processes you want
+		while (process != processors) {
+			int pid = fork();
+			
+			if (pid > 0) {
+				processIDS[lines[process]->end] = pid;  //create map from line number to pid so you can append files in correct order later
+				process++;
+			}else if (pid == 0){
+				driver(lines[process]->start, lines[process]->end, filename + toString(getpid()) + ".temp", cutoff);
+				exit(0);
+			}else { cout << "unable to spawn the necessary processes." << endl; exit(0); }
+		}
+	
+		//force parent to wait until all the processes are done
+		for (map<int, int>::iterator it = processIDS.begin(); it != processIDS.end(); it++) { 
+			int temp = it->second;
+			wait(&temp);
+		}
+#endif
+	}
+	catch(exception& e) {
+		cout << "Standard Error: " << e.what() << " has occurred in the DistanceCommand class Function createProcesses. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
+		exit(1);
+	}
+	catch(...) {
+		cout << "An unknown error has occurred in the DistanceCommand class function createProcesses. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
+		exit(1);
+	}	
+}
 
 /**************************************************************************************************/
 /////// need to fix to work with calcs and sequencedb
-int DistanceCommand::driver(Dist* distCalculator, SequenceDB* align, int startLine, int endLine, string dFileName, float cutoff){
+int DistanceCommand::driver(int startLine, int endLine, string dFileName, float cutoff){
 	try {
+
 		int startTime = time(NULL);
-	
-		ofstream distFile(dFileName.c_str(), ios::trunc);
-		distFile.setf(ios::fixed, ios::showpoint);
-		distFile << setprecision(4);
-	
-		for(int i=startLine;i<endLine;i++){
 		
+		//column file
+		ofstream outFile(dFileName.c_str(), ios::trunc);
+		outFile.setf(ios::fixed, ios::showpoint);
+		outFile << setprecision(4);
+		
+		if(isTrue(phylip) && startLine == 0){	outFile << alignDB.getNumSeqs() << endl;	}
+		for(int i=startLine;i<endLine;i++){
+			if(isTrue(phylip))	{	outFile << alignDB.get(i).getName() << '\t';	}
 			for(int j=0;j<i;j++){
-				distCalculator->calcDist(align->get(i), align->get(j));
+				distCalculator->calcDist(alignDB.get(i), alignDB.get(j));
 				double dist = distCalculator->getDist();
-
+				
 				if(dist <= cutoff){
-					distFile << align->get(i).getName() << ' ' << align->get(j).getName() << ' ' << dist << endl;
+					if (!isTrue(phylip)) { outFile << alignDB.get(i).getName() << ' ' << alignDB.get(j).getName() << ' ' << dist << endl; }
 				}
-			
+				if (isTrue(phylip)) {  outFile << dist << '\t'; }
+				
 			}
+			
+			if (isTrue(phylip) == true) { outFile << endl; }
+			
 			if(i % 100 == 0){
 				cout << i << '\t' << time(NULL) - startTime << endl;
 			}
-		
+			
 		}
 		cout << endLine-1 << '\t' << time(NULL) - startTime << endl;
-	
+		
+		outFile.close();
+		
 		return 1;
 	}
 	catch(exception& e) {
@@ -227,7 +297,7 @@ int DistanceCommand::driver(Dist* distCalculator, SequenceDB* align, int startLi
 		cout << "An unknown error has occurred in the DistanceCommand class function driver. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
 		exit(1);
 	}	
-
+	
 }
 
 /**************************************************************************************************/
@@ -235,7 +305,7 @@ void DistanceCommand::appendFiles(string temp, string filename) {
 	try{
 		ofstream output;
 		ifstream input;
-		
+	
 		//open output file in append mode
 		openOutputFileAppend(filename, output);
 		
@@ -250,7 +320,7 @@ void DistanceCommand::appendFiles(string temp, string filename) {
 				output << line << endl;   // Appending back newline char 
 			}
 		}	
-			
+		
 		input.close();
 		output.close();
 	}