]> git.donarmstrong.com Git - mothur.git/blobdiff - distancecommand.cpp
working on megastats
[mothur.git] / distancecommand.cpp
index 351a3f4f758a28bb78f8277e520ae74092e5b740..871e07b65510d33214d78c9efee3b3448d10d405 100644 (file)
@@ -26,7 +26,8 @@ DistanceCommand::DistanceCommand(string option) {
                
                else {
                        //valid paramters for this command
-                       string Array[] =  {"fasta","oldfasta","column", "output", "calc", "countends", "cutoff", "processors", "outputdir","inputdir"};
+                       string Array[] =  {"fasta","oldfasta","column", "output", "calc", "countends", "cutoff", "processors", "outputdir","inputdir","compress"};
+
                        vector<string> myArray (Array, Array+(sizeof(Array)/sizeof(string)));
                        
                        OptionParser parser(option);
@@ -114,6 +115,9 @@ DistanceCommand::DistanceCommand(string option) {
                        temp = validParameter.validFile(parameters, "processors", false);       if(temp == "not found"){        temp = "1"; }
                        convert(temp, processors); 
                        
+                       temp = validParameter.validFile(parameters, "compress", false);         if(temp == "not found"){  temp = "F"; }
+                       convert(temp, compress);
+
                        output = validParameter.validFile(parameters, "output", false);         if(output == "not found"){      output = "column"; }
                        
                        if (((column != "") && (oldfastafile == "")) || ((column == "") && (oldfastafile != ""))) { m->mothurOut("If you provide column or oldfasta, you must provide both."); m->mothurOutEndLine(); abort=true; }
@@ -166,7 +170,7 @@ DistanceCommand::~DistanceCommand(){
 void DistanceCommand::help(){
        try {
                m->mothurOut("The dist.seqs command reads a file containing sequences and creates a distance file.\n");
-               m->mothurOut("The dist.seqs command parameters are fasta, oldfasta, column, calc, countends, output, cutoff and processors.  \n");
+               m->mothurOut("The dist.seqs command parameters are fasta, oldfasta, column, calc, countends, output, compress, cutoff and processors.  \n");
                m->mothurOut("The fasta parameter is required.\n");
                m->mothurOut("The oldfasta and column parameters allow you to append the distances calculated to the column file.\n");
                m->mothurOut("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");
@@ -174,6 +178,7 @@ void DistanceCommand::help(){
                m->mothurOut("The cutoff parameter allows you to specify maximum distance to keep. The default is 1.0.\n");
                m->mothurOut("The output parameter allows you to specify format of your distance matrix. Options are column, lt, and square. The default is column.\n");
                m->mothurOut("The processors parameter allows you to specify number of processors to use.  The default is 1.\n");
+               m->mothurOut("The compress parameter allows you to indicate that you want the resulting distance file compressed.  The default is false.\n");
                m->mothurOut("The dist.seqs command should be in the following format: \n");
                m->mothurOut("dist.seqs(fasta=yourFastaFile, calc=yourCalc, countends=yourEnds, cutoff= yourCutOff, processors=yourProcessors) \n");
                m->mothurOut("Example dist.seqs(fasta=amazon.fasta, calc=eachgap, countends=F, cutoff= 2.0, processors=3).\n");
@@ -237,10 +242,15 @@ int DistanceCommand::execute(){
                MPI_Comm_rank(MPI_COMM_WORLD, &pid); //find out who we are
                
                //each process gets where it should start and stop in the file
-               start = int (sqrt(float(pid)/float(processors)) * numSeqs);
-               end = int (sqrt(float(pid+1)/float(processors)) * numSeqs);
+               if (output != "square") {
+                       start = int (sqrt(float(pid)/float(processors)) * numSeqs);
+                       end = int (sqrt(float(pid+1)/float(processors)) * numSeqs);
+               }else{
+                       start = int ((float(pid)/float(processors)) * numSeqs);
+                       end = int ((float(pid+1)/float(processors)) * numSeqs);
+               }
                
-               if (output != "lt") {
+               if (output == "column") {
                        MPI_File outMPI;
                        int amode=MPI_MODE_CREATE|MPI_MODE_WRONLY; 
 
@@ -257,7 +267,8 @@ int DistanceCommand::execute(){
                        
                                //do your part
                                string outputMyPart;
-                               driverMPI(start, end, outMPI, cutoff);
+                               
+                               driverMPI(start, end, outMPI, cutoff); 
                                
                                if (m->control_pressed) { MPI_File_close(&outMPI);  delete distCalculator;  return 0; }
                        
@@ -270,7 +281,7 @@ int DistanceCommand::execute(){
                                }
                        }else { //you are a child process
                                //do your part
-                               driverMPI(start, end, outMPI, cutoff);
+                               driverMPI(start, end, outMPI, cutoff); 
                                
                                if (m->control_pressed) { MPI_File_close(&outMPI);  delete distCalculator;  return 0; }
                        
@@ -287,8 +298,10 @@ int DistanceCommand::execute(){
                        
                                //do your part
                                string outputMyPart;
-                               long mySize;
-                               driverMPI(start, end, outputFile, mySize);
+                               unsigned long int mySize;
+                               
+                               if (output != "square"){ driverMPI(start, end, outputFile, mySize); }
+                               else { driverMPI(start, end, outputFile, mySize, output); }
        
                                if (m->control_pressed) {  delete distCalculator;  return 0; }
                                
@@ -307,7 +320,7 @@ int DistanceCommand::execute(){
 
                                //wait on chidren
                                for(int b = 1; b < processors; b++) { 
-                                       long fileSize;
+                                       unsigned long int fileSize;
                                        
                                        if (m->control_pressed) { MPI_File_close(&outMPI);  delete distCalculator;  return 0; }
                                        
@@ -335,8 +348,9 @@ int DistanceCommand::execute(){
                                MPI_File_close(&outMPI);
                        }else { //you are a child process
                                //do your part
-                               long size;
-                               driverMPI(start, end, (outputFile + toString(pid) + ".temp"), size);
+                               unsigned long int size;
+                               if (output != "square"){ driverMPI(start, end, (outputFile + toString(pid) + ".temp"), size); }
+                               else { driverMPI(start, end, (outputFile + toString(pid) + ".temp"), size, output); }
                                
                                if (m->control_pressed) { delete distCalculator;  return 0; }
                        
@@ -350,13 +364,19 @@ int DistanceCommand::execute(){
        #if defined (__APPLE__) || (__MACH__) || (linux) || (__linux)
                //if you don't need to fork anything
                if(processors == 1){
-                       driver(0, numSeqs, outputFile, cutoff);
+                       if (output != "square") {  driver(0, numSeqs, outputFile, cutoff); }
+                       else { driver(0, numSeqs, outputFile, "square");  }
                }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);
+                               if (output != "square") {
+                                       lines[i]->start = int (sqrt(float(i)/float(processors)) * numSeqs);
+                                       lines[i]->end = int (sqrt(float(i+1)/float(processors)) * numSeqs);
+                               }else{
+                                       lines[i]->start = int ((float(i)/float(processors)) * numSeqs);
+                                       lines[i]->end = int ((float(i+1)/float(processors)) * numSeqs);
+                               }
                        }
 
                        createProcesses(outputFile); 
@@ -372,8 +392,9 @@ int DistanceCommand::execute(){
                        }
                }
        #else
-               ifstream inFASTA;
-               driver(0, numSeqs, outputFile, cutoff);
+               //ifstream inFASTA;
+               if (output != "square") {  driver(0, numSeqs, outputFile, cutoff); }
+               else { driver(0, numSeqs, outputFile, "square");  }
        #endif
        
 #endif
@@ -385,7 +406,7 @@ int DistanceCommand::execute(){
                        if (pid == 0) { //only one process should output to screen
                #endif
                
-               if (output == "square") {  convertMatrix(outputFile); }
+               //if (output == "square") {  convertMatrix(outputFile); }
                
                ifstream fileHandle;
                fileHandle.open(outputFile.c_str());
@@ -406,6 +427,13 @@ int DistanceCommand::execute(){
                                remove(outputFile.c_str());
                                outputFile = column;
                        }
+                       
+                       if (outputDir != "") { 
+                               string newOutputName = outputDir + m->getSimpleName(outputFile);
+                               rename(outputFile.c_str(), newOutputName.c_str());
+                               remove(outputFile.c_str());
+                               outputFile = newOutputName;
+                       }
                }
 
                
@@ -422,6 +450,14 @@ int DistanceCommand::execute(){
                m->mothurOut(outputFile); m->mothurOutEndLine();
                m->mothurOutEndLine();
                m->mothurOut("It took " + toString(time(NULL) - startTime) + " to calculate the distances for " + toString(numSeqs) + " sequences."); m->mothurOutEndLine();
+
+
+               if (m->isTrue(compress)) {
+                       m->mothurOut("Compressing..."); m->mothurOutEndLine();
+                       m->mothurOut("(Replacing " + outputFile + " with " + outputFile + ".gz)"); m->mothurOutEndLine();
+                       system(("gzip -v " + outputFile).c_str());
+               }
+
                return 0;
                
        }
@@ -445,7 +481,8 @@ void DistanceCommand::createProcesses(string filename) {
                                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);
+                               if (output != "square") {  driver(lines[process]->start, lines[process]->end, filename + toString(getpid()) + ".temp", cutoff); }
+                               else { driver(lines[process]->start, lines[process]->end, filename + toString(getpid()) + ".temp", "square"); }
                                exit(0);
                        }else { m->mothurOut("unable to spawn the necessary processes."); m->mothurOutEndLine(); exit(0); }
                }
@@ -500,12 +537,6 @@ int DistanceCommand::driver(int startLine, int endLine, string dFileName, float
                                        if (output == "column") { outFile << alignDB.get(i).getName() << ' ' << alignDB.get(j).getName() << ' ' << dist << endl; }
                                }
                                if (output == "lt") {  outFile << dist << '\t'; }
-                               
-                               if (output == "square") { //make a square column you can convert to square phylip
-                                       outFile << alignDB.get(i).getName() << '\t' << alignDB.get(j).getName() << '\t' << dist << endl;
-                                       outFile << alignDB.get(j).getName() << '\t' << alignDB.get(i).getName() << '\t' << dist << endl;
-                               }
-
                        }
                        
                        if (output == "lt") { outFile << endl; }
@@ -526,6 +557,56 @@ int DistanceCommand::driver(int startLine, int endLine, string dFileName, float
                exit(1);
        }
 }
+/**************************************************************************************************/
+/////// need to fix to work with calcs and sequencedb
+int DistanceCommand::driver(int startLine, int endLine, string dFileName, string square){
+       try {
+
+               int startTime = time(NULL);
+               
+               //column file
+               ofstream outFile(dFileName.c_str(), ios::trunc);
+               outFile.setf(ios::fixed, ios::showpoint);
+               outFile << setprecision(4);
+               
+               if(startLine == 0){     outFile << alignDB.getNumSeqs() << endl;        }
+               
+               for(int i=startLine;i<endLine;i++){
+                               
+                       string name = alignDB.get(i).getName();
+                       //pad with spaces to make compatible
+                       if (name.length() < 10) { while (name.length() < 10) {  name += " ";  } }
+                               
+                       outFile << name << '\t';        
+                       
+                       for(int j=0;j<alignDB.getNumSeqs();j++){
+                               
+                               if (m->control_pressed) { outFile.close(); return 0;  }
+                               
+                               distCalculator->calcDist(alignDB.get(i), alignDB.get(j));
+                               double dist = distCalculator->getDist();
+                               
+                               outFile << dist << '\t'; 
+                       }
+                       
+                       outFile << endl; 
+                       
+                       if(i % 100 == 0){
+                               m->mothurOut(toString(i) + "\t" + toString(time(NULL) - startTime)); m->mothurOutEndLine();
+                       }
+                       
+               }
+               m->mothurOut(toString(endLine-1) + "\t" + toString(time(NULL) - startTime)); m->mothurOutEndLine();
+               
+               outFile.close();
+               
+               return 1;
+       }
+       catch(exception& e) {
+               m->errorOut(e, "DistanceCommand", "driver");
+               exit(1);
+       }
+}
 #ifdef USE_MPI
 /**************************************************************************************************/
 /////// need to fix to work with calcs and sequencedb
@@ -550,12 +631,7 @@ int DistanceCommand::driverMPI(int startLine, int endLine, MPI_File& outMPI, flo
                                double dist = distCalculator->getDist();
                                
                                if(dist <= cutoff){
-                                       if (output == "column") { outputString += (alignDB.get(i).getName() + ' ' + alignDB.get(j).getName() + ' ' + toString(dist) + '\n'); }
-                               }
-                               
-                               if (output == "square") { //make a square column you can convert to square phylip
-                                       outputString += (alignDB.get(i).getName() + ' ' + alignDB.get(j).getName() + ' ' + toString(dist) + '\n');
-                                       outputString += (alignDB.get(j).getName() + ' ' + alignDB.get(i).getName() + ' ' + toString(dist) + '\n');
+                                        outputString += (alignDB.get(i).getName() + ' ' + alignDB.get(j).getName() + ' ' + toString(dist) + '\n'); 
                                }
                        }
                        
@@ -588,7 +664,7 @@ int DistanceCommand::driverMPI(int startLine, int endLine, MPI_File& outMPI, flo
 }
 /**************************************************************************************************/
 /////// need to fix to work with calcs and sequencedb
-int DistanceCommand::driverMPI(int startLine, int endLine, string file, long& size){
+int DistanceCommand::driverMPI(int startLine, int endLine, string file, unsigned long int& size){
        try {
                MPI_Status status;
                
@@ -609,16 +685,16 @@ int DistanceCommand::driverMPI(int startLine, int endLine, string file, long& si
                string outputString = "";
                size = 0;
                
-               if((output == "lt") && startLine == 0){ outputString += toString(alignDB.getNumSeqs()) + "\n";  }
+               if(startLine == 0){     outputString += toString(alignDB.getNumSeqs()) + "\n";  }
                
                for(int i=startLine;i<endLine;i++){
-                       if(output == "lt")      {       
-                               string name = alignDB.get(i).getName();
-                               if (name.length() < 10) { //pad with spaces to make compatible
-                                       while (name.length() < 10) {  name += " ";  }
-                               }
-                               outputString += name + "\t";    
+                               
+                       string name = alignDB.get(i).getName();
+                       if (name.length() < 10) { //pad with spaces to make compatible
+                               while (name.length() < 10) {  name += " ";  }
                        }
+                       outputString += name + "\t";    
+                       
                        for(int j=0;j<i;j++){
                                
                                if (m->control_pressed) {  return 0;  }
@@ -626,10 +702,10 @@ int DistanceCommand::driverMPI(int startLine, int endLine, string file, long& si
                                distCalculator->calcDist(alignDB.get(i), alignDB.get(j));
                                double dist = distCalculator->getDist();
                                
-                               if (output == "lt") {  outputString += toString(dist) + "\t"; }
+                               outputString += toString(dist) + "\t"; 
                        }
                        
-                       if (output == "lt") { outputString += "\n"; }
+                       outputString += "\n"; 
 
                
                        if(i % 100 == 0){
@@ -660,8 +736,82 @@ int DistanceCommand::driverMPI(int startLine, int endLine, string file, long& si
                exit(1);
        }
 }
-#endif
 /**************************************************************************************************/
+/////// need to fix to work with calcs and sequencedb
+int DistanceCommand::driverMPI(int startLine, int endLine, string file, unsigned long int& size, string square){
+       try {
+               MPI_Status status;
+               
+               MPI_File outMPI;
+               int amode=MPI_MODE_CREATE|MPI_MODE_WRONLY; 
+
+               //char* filename = new char[file.length()];
+               //memcpy(filename, file.c_str(), file.length());
+               
+               char filename[1024];
+               strcpy(filename, file.c_str());
+
+               MPI_File_open(MPI_COMM_SELF, filename, amode, MPI_INFO_NULL, &outMPI);
+               //delete filename;
+
+               int startTime = time(NULL);
+               
+               string outputString = "";
+               size = 0;
+               
+               if(startLine == 0){     outputString += toString(alignDB.getNumSeqs()) + "\n";  }
+               
+               for(int i=startLine;i<endLine;i++){
+                               
+                       string name = alignDB.get(i).getName();
+                       if (name.length() < 10) { //pad with spaces to make compatible
+                               while (name.length() < 10) {  name += " ";  }
+                       }
+                       outputString += name + "\t";    
+                       
+                       for(int j=0;j<alignDB.getNumSeqs();j++){
+                               
+                               if (m->control_pressed) {  return 0;  }
+                               
+                               distCalculator->calcDist(alignDB.get(i), alignDB.get(j));
+                               double dist = distCalculator->getDist();
+                               
+                               outputString += toString(dist) + "\t"; 
+                       }
+                       
+                       outputString += "\n"; 
+
+               
+                       if(i % 100 == 0){
+                               //m->mothurOut(toString(i) + "\t" + toString(time(NULL) - startTime)); m->mothurOutEndLine();
+                               cout << i << '\t' << (time(NULL) - startTime) << endl;
+                       }
+                       
+                       
+                       //send results to parent
+                       int length = outputString.length();
+                       char* buf = new char[length];
+                       memcpy(buf, outputString.c_str(), length);
+                       
+                       MPI_File_write(outMPI, buf, length, MPI_CHAR, &status);
+                       size += outputString.length();
+                       outputString = "";
+                       delete buf;
+               }
+               
+               //m->mothurOut(toString(endLine-1) + "\t" + toString(time(NULL) - startTime)); m->mothurOutEndLine();
+               cout << (endLine-1) << '\t' << (time(NULL) - startTime) << endl;
+               MPI_File_close(&outMPI);
+               
+               return 1;
+       }
+       catch(exception& e) {
+               m->errorOut(e, "DistanceCommand", "driverMPI");
+               exit(1);
+       }
+}
+#endif
+/**************************************************************************************************
 int DistanceCommand::convertMatrix(string outputFile) {
        try{
 
@@ -749,7 +899,7 @@ int DistanceCommand::convertMatrix(string outputFile) {
                exit(1);
        }
 }
-/**************************************************************************************************/
+/**************************************************************************************************
 int DistanceCommand::convertToLowerTriangle(string outputFile) {
        try{
 
@@ -933,30 +1083,8 @@ bool DistanceCommand::sanityCheck() {
                exit(1);
        }
 }
-
-/**************************************************************************************************
-void DistanceCommand::m->appendFiles(string temp, string filename) {
-       try{
-               ofstream output;
-               ifstream input;
-       
-               //open output file in append mode
-               m->openOutputFileAppend(filename, output);
-               m->openInputFile(temp, input);
-               
-               while(char c = input.get()){
-                       if(input.eof())         {       break;                  }
-                       else                            {       output << c;    }
-               }
-               
-               input.close();
-               output.close();
-       }
-       catch(exception& e) {
-               m->errorOut(e, "DistanceCommand", "m->appendFiles");
-               exit(1);
-       }
-}
 /**************************************************************************************************/
 
 
+
+