X-Git-Url: https://git.donarmstrong.com/?a=blobdiff_plain;f=distancecommand.cpp;h=7890cb9303cd39e38967a7694203fac8d2d588ab;hb=0470f6d037aacb3563c3f7010708120a4a67d4e6;hp=0714ab3623181a86d5e0485b5a4e6613144a6b66;hpb=d5d2761f88b41f1006d0b700e0ab51e2ce48b875;p=mothur.git diff --git a/distancecommand.cpp b/distancecommand.cpp index 0714ab3..7890cb9 100644 --- a/distancecommand.cpp +++ b/distancecommand.cpp @@ -16,44 +16,93 @@ //********************************************************************************************************************** -DistanceCommand::DistanceCommand(){ +DistanceCommand::DistanceCommand(string option){ try { globaldata = GlobalData::getInstance(); - validCalculator = new ValidCalculators(); - ends = globaldata->getEnds(); - convert(globaldata->getProcessors(), processors); - convert(globaldata->getCutOff(), cutoff); + abort = false; + Estimators.clear(); - int i; - if (ends != "T") { - for (i=0; iEstimators.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") { + //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 myArray (Array, Array+(sizeof(Array)/sizeof(string))); + + parser = new OptionParser(); + parser->parse(option, parameters); delete parser; + + ValidParameters* validParameter = new ValidParameters(); + + //check to make sure all parameters are valid for command + for (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 { + globaldata->setFastaFile(fastafile); + openInputFile(fastafile, in); + } + + //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"; } + + delete validParameter; + + validCalculator = new ValidCalculators(); + + int i; + if (isTrue(countends) == true) { + for (i=0; iisValidCalculator("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 (i=0; iEstimators.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(); + }else { + for (i=0; iisValidCalculator("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(); + } } } } + + delete validCalculator; } - //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"; @@ -66,118 +115,101 @@ DistanceCommand::DistanceCommand(){ } //********************************************************************************************************************** +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(); - + //reads fasta file and fills sequenceDB + seqDB = new SequenceDB(in); + int numSeqs = seqDB->getNumSeqs(); + cutoff += 0.005; - string distFile = getRootName(globaldata->getFastaFile()) + "dist"; + string outputFile; - remove(distFile.c_str()); + //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 + openOutputFile(outputFile, outFile); + outFile << numSeqs << endl; + outFile.close(); + }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, phylipFile, cutoff); + //# endif + +#if defined (__APPLE__) || (__MACH__) || (linux) || (__linux) + //if you don't need to fork anything 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()); + driver(distCalculator, seqDB, 0, numSeqs, outputFile + ".temp", cutoff); + appendFiles((outputFile + ".temp"), outputFile); + remove((outputFile + ".temp").c_str()); + }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); } - 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); - } - 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); - //system(("cat " + distFile + "tempa" + " >> " + distFile).c_str()); - appendFiles(distFile+"tempa", distFile); - //system(("rm " + distFile + "tempa").c_str()); - remove((distFile + "tempa").c_str()); - } - else{ - driver(distCalculator, seqDB, numSeqs / 2, (numSeqs/sqrt(2)), distFile + "tempb", cutoff); - //system(("cat " + distFile + "tempb" + " >> " + distFile).c_str()); - appendFiles(distFile+"tempb", distFile); - //system(("rm " + distFile + "tempb").c_str()); - 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); - //system(("cat " + distFile + "tempc" + " >> " + distFile).c_str()); - appendFiles(distFile+"tempc", distFile); - //system(("rm " + distFile + "tempc").c_str()); - remove((distFile + "tempc").c_str()); - } - else{ - driver(distCalculator, seqDB, (sqrt(3) * numSeqs / 2), numSeqs, distFile + "tempd", cutoff); - //system(("cat " + distFile + "tempd" + " >> " + distFile).c_str()); - appendFiles(distFile+"tempd", distFile); - //system(("rm " + distFile + "tempd").c_str()); - remove((distFile + "tempd").c_str()); - } - wait(NULL); + cout << lines[0]->start << '\t' << lines[0]->end << endl; + cout << lines[1]->start << '\t' << lines[1]->end << endl; + + createProcesses(outputFile); + + //append and remove temp files + for (it = processIDS.begin(); it != processIDS.end(); it++) { + appendFiles((outputFile + toString(it->second) + ".temp"), outputFile); + remove((outputFile + toString(it->second) + ".temp").c_str()); } - wait(NULL); } - wait(NULL); - +#else + driver(distCalculator, seqDB, 0, numSeqs, outputFile + ".temp", cutoff); + appendFiles((outputFile + ".temp"), outputFile); + remove((outputFile + ".temp").c_str()); +#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,36 +220,79 @@ int DistanceCommand::execute(){ exit(1); } } +/**************************************************************************************************/ +void DistanceCommand::createProcesses(string filename) { + try { + 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(distCalculator, seqDB, 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 (it = processIDS.begin(); it != processIDS.end(); it++) { + int temp = it->second; + wait(&temp); + } + + } + 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){ 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;icalcDist(align->get(i), align->get(j)); + //column file + ofstream outFile(dFileName.c_str(), ios::trunc); + outFile.setf(ios::fixed, ios::showpoint); + outFile << setprecision(4); + + for(int i=startLine;icalcDist(*(align->get(i)), *(align->get(j))); double dist = distCalculator->getDist(); - + if(dist <= cutoff){ - distFile << align->get(i).getName() << ' ' << align->get(j).getName() << ' ' << dist << endl; -//cout << align->get(i).getName() << ' ' << align->get(j).getName() << ' ' << dist << endl; + if (isTrue(phylip) != true) { outFile << align->get(i)->getName() << ' ' << align->get(j)->getName() << ' ' << dist << endl; } } - + if (isTrue(phylip) == true) { 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; - + + //philFile.close(); + //distFile.close(); + return 1; } catch(exception& e) { @@ -228,7 +303,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); } - + } /**************************************************************************************************/ @@ -236,7 +311,7 @@ void DistanceCommand::appendFiles(string temp, string filename) { try{ ofstream output; ifstream input; - + //open output file in append mode openOutputFileAppend(filename, output); @@ -251,7 +326,7 @@ void DistanceCommand::appendFiles(string temp, string filename) { output << line << endl; // Appending back newline char } } - + input.close(); output.close(); } @@ -264,3 +339,4 @@ void DistanceCommand::appendFiles(string temp, string filename) { exit(1); } } +/**************************************************************************************************/ \ No newline at end of file