]> git.donarmstrong.com Git - mothur.git/blobdiff - distancecommand.cpp
paralellized screen.seqs and added mpi code to it. fixed bug with all mpi commands...
[mothur.git] / distancecommand.cpp
index 13f56eb90a422d43e821a9fa378f582691468cb8..ec32ffeb9e2e28b7398f78c5dd5fd858de2326a0 100644 (file)
@@ -162,11 +162,13 @@ int DistanceCommand::execute(){
                
                if (abort == true) { return 0; }
                
+               int startTime = time(NULL);
+               
                int numSeqs = alignDB.getNumSeqs();
                cutoff += 0.005;
                
                string outputFile;
-               
+                               
                if (output == "lt") { //does the user want lower triangle phylip formatted file 
                        outputFile = outputDir + getRootName(getSimpleName(fastafile)) + "phylip.dist";
                        remove(outputFile.c_str());
@@ -179,8 +181,128 @@ int DistanceCommand::execute(){
                        outputFile = outputDir + getRootName(getSimpleName(fastafile)) + "square.dist";
                        remove(outputFile.c_str());
                }
+               
+
+#ifdef USE_MPI
+               
+               int pid, start, end; 
+               int tag = 2001;
                                
-#if defined (__APPLE__) || (__MACH__) || (linux) || (__linux)
+               MPI_Status status; 
+               MPI_Comm_size(MPI_COMM_WORLD, &processors); //set processors to the number of mpi processes running
+               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 != "lt") {
+                       MPI_File outMPI;
+                       int amode=MPI_MODE_CREATE|MPI_MODE_WRONLY; 
+
+                       //char* filename = new char[outputFile.length()];
+                       //memcpy(filename, outputFile.c_str(), outputFile.length());
+                       
+                       char filename[1024];
+                       strcpy(filename, outputFile.c_str());
+                       
+                       MPI_File_open(MPI_COMM_WORLD, filename, amode, MPI_INFO_NULL, &outMPI);
+                       //delete filename;
+
+                       if (pid == 0) { //you are the root process 
+                       
+                               //do your part
+                               string outputMyPart;
+                               driverMPI(start, end, outMPI, cutoff);
+                               
+                               if (m->control_pressed) { MPI_File_close(&outMPI);  delete distCalculator;  return 0; }
+                       
+                               //wait on chidren
+                               for(int i = 1; i < processors; i++) { 
+                                       if (m->control_pressed) { MPI_File_close(&outMPI);  delete distCalculator;  return 0; }
+                                       
+                                       char buf[4];
+                                       MPI_Recv(buf, 4, MPI_CHAR, i, tag, MPI_COMM_WORLD, &status); 
+                               }
+                       }else { //you are a child process
+                               //do your part
+                               driverMPI(start, end, outMPI, cutoff);
+                               
+                               if (m->control_pressed) { MPI_File_close(&outMPI);  delete distCalculator;  return 0; }
+                       
+                               char buf[4];
+                               strcpy(buf, "done"); 
+                               //tell parent you are done.
+                               MPI_Send(buf, 4, MPI_CHAR, 0, tag, MPI_COMM_WORLD);
+                       }
+                       
+                       MPI_File_close(&outMPI);
+                       
+               }else { //lower triangle format
+                       if (pid == 0) { //you are the root process 
+                       
+                               //do your part
+                               string outputMyPart;
+                               long mySize;
+                               driverMPI(start, end, outputFile, mySize);
+       
+                               if (m->control_pressed) {  delete distCalculator;  return 0; }
+                               
+                               int amode=MPI_MODE_APPEND|MPI_MODE_WRONLY|MPI_MODE_CREATE; //
+                               MPI_File outMPI;
+                               MPI_File inMPI;
+
+                               //char* filename = new char[outputFile.length()];
+                               //memcpy(filename, outputFile.c_str(), outputFile.length());
+                               
+                               char filename[1024];
+                               strcpy(filename, outputFile.c_str());
+
+                               MPI_File_open(MPI_COMM_SELF, filename, amode, MPI_INFO_NULL, &outMPI);
+                               //delete filename;
+
+                               //wait on chidren
+                               for(int b = 1; b < processors; b++) { 
+                                       long fileSize;
+                                       
+                                       if (m->control_pressed) { MPI_File_close(&outMPI);  delete distCalculator;  return 0; }
+                                       
+                                       MPI_Recv(&fileSize, 1, MPI_LONG, b, tag, MPI_COMM_WORLD, &status); 
+                                       
+                                       string outTemp = outputFile + toString(b) + ".temp";
+
+                                       char* buf = new char[outTemp.length()];
+                                       memcpy(buf, outTemp.c_str(), outTemp.length());
+                                       
+                                       MPI_File_open(MPI_COMM_SELF, buf, MPI_MODE_DELETE_ON_CLOSE|MPI_MODE_RDONLY, MPI_INFO_NULL, &inMPI);
+                                       delete buf;
+
+                                       int count = 0;
+                                       while (count < fileSize) { 
+                                               char buf2[1];
+                                               MPI_File_read(inMPI, buf2, 1, MPI_CHAR, &status);
+                                               MPI_File_write(outMPI, buf2, 1, MPI_CHAR, &status);
+                                               count += 1;
+                                       }
+                                       
+                                       MPI_File_close(&inMPI); //deleted on close
+                               }
+                               
+                               MPI_File_close(&outMPI);
+                       }else { //you are a child process
+                               //do your part
+                               long size;
+                               driverMPI(start, end, (outputFile + toString(pid) + ".temp"), size);
+                               
+                               if (m->control_pressed) { delete distCalculator;  return 0; }
+                       
+                               //tell parent you are done.
+                               MPI_Send(&size, 1, MPI_LONG, 0, tag, MPI_COMM_WORLD);
+                       }
+               }
+#else          
+                               
+       #if defined (__APPLE__) || (__MACH__) || (linux) || (__linux)
                //if you don't need to fork anything
                if(processors == 1){
                        driver(0, numSeqs, outputFile, cutoff);
@@ -204,20 +326,43 @@ int DistanceCommand::execute(){
                                remove((outputFile + toString(it->second) + ".temp").c_str());
                        }
                }
-#else
+       #else
                ifstream inFASTA;
                driver(0, numSeqs, outputFile, cutoff);
+       #endif
+       
 #endif
+               if (m->control_pressed) { delete distCalculator; remove(outputFile.c_str()); return 0; }
+               
+               #ifdef USE_MPI
+                       MPI_Comm_rank(MPI_COMM_WORLD, &pid); 
+                                       
+                       if (pid == 0) { //only one process should output to screen
+               #endif
                
                if (output == "square") {  convertMatrix(outputFile); }
                
+               ifstream fileHandle;
+               fileHandle.open(outputFile.c_str());
+               if(fileHandle) {
+                       gobble(fileHandle);
+                       if (fileHandle.eof()) { m->mothurOut(outputFile + " is blank. This can result if there are no distances below your cutoff.");  m->mothurOutEndLine(); }
+               }
+
+               
+               #ifdef USE_MPI
+                       }
+               #endif
+               
+               if (m->control_pressed) { delete distCalculator; remove(outputFile.c_str()); return 0; }
+               
                delete distCalculator;
                
                m->mothurOutEndLine();
                m->mothurOut("Output File Name: "); m->mothurOutEndLine();
                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();
                return 0;
                
        }
@@ -282,6 +427,9 @@ int DistanceCommand::driver(int startLine, int endLine, string dFileName, float
                                outFile << name << '\t';        
                        }
                        for(int j=0;j<i;j++){
+                               
+                               if (m->control_pressed) { outFile.close(); return 0;  }
+                               
                                distCalculator->calcDist(alignDB.get(i), alignDB.get(j));
                                double dist = distCalculator->getDist();
                                
@@ -315,8 +463,139 @@ int DistanceCommand::driver(int startLine, int endLine, string dFileName, float
                exit(1);
        }
 }
+#ifdef USE_MPI
+/**************************************************************************************************/
+/////// need to fix to work with calcs and sequencedb
+int DistanceCommand::driverMPI(int startLine, int endLine, MPI_File& outMPI, float cutoff){
+       try {
+               MPI_Status status;
+               int startTime = time(NULL);
+               
+               string outputString = "";
+               
+               for(int i=startLine;i<endLine;i++){
+       
+                       for(int j=0;j<i;j++){
+                               
+                               if (m->control_pressed) {  return 0;  }
+                               
+                               distCalculator->calcDist(alignDB.get(i), alignDB.get(j));
+                               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');
+                               }
+                       }
+                       
+                       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_shared(outMPI, buf, length, MPI_CHAR, &status);
+                       outputString = "";
+                       delete buf;
+                       
+               }
+               
+               //m->mothurOut(toString(endLine-1) + "\t" + toString(time(NULL) - startTime)); m->mothurOutEndLine();
+               cout << (endLine-1) << '\t' << (time(NULL) - startTime) << endl;                
+               return 1;
+       }
+       catch(exception& e) {
+               m->errorOut(e, "DistanceCommand", "driverMPI");
+               exit(1);
+       }
+}
+/**************************************************************************************************/
+/////// need to fix to work with calcs and sequencedb
+int DistanceCommand::driverMPI(int startLine, int endLine, string file, long& size){
+       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((output == "lt") && 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";    
+                       }
+                       for(int j=0;j<i;j++){
+                               
+                               if (m->control_pressed) {  return 0;  }
+                               
+                               distCalculator->calcDist(alignDB.get(i), alignDB.get(j));
+                               double dist = distCalculator->getDist();
+                               
+                               if (output == "lt") {  outputString += toString(dist) + "\t"; }
+                       }
+                       
+                       if (output == "lt") { 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
 /**************************************************************************************************/
-void DistanceCommand::convertMatrix(string outputFile) {
+int DistanceCommand::convertMatrix(string outputFile) {
        try{
 
                //sort file by first column so the distances for each row are together
@@ -358,6 +637,8 @@ void DistanceCommand::convertMatrix(string outputFile) {
                //openInputFile(outfile, in);
                
                while(!in.eof()) {
+                       if (m->control_pressed) { in.close(); remove(outfile.c_str()); out.close(); return 0; }
+                       
                        in >> first >> second >> dist; gobble(in);
                                
                        if (first != currentRow) {
@@ -393,12 +674,108 @@ void DistanceCommand::convertMatrix(string outputFile) {
                
                remove(outfile.c_str());
                
+               return 1;
+               
        }
        catch(exception& e) {
                m->errorOut(e, "DistanceCommand", "convertMatrix");
                exit(1);
        }
 }
+/**************************************************************************************************/
+int DistanceCommand::convertToLowerTriangle(string outputFile) {
+       try{
+
+               //sort file by first column so the distances for each row are together
+               string outfile = getRootName(outputFile) + "sorted.dist.temp";
+               
+               //use the unix sort 
+               #if defined (__APPLE__) || (__MACH__) || (linux) || (__linux)
+                       string command = "sort -n " + outputFile + " -o " + outfile;
+                       system(command.c_str());
+               #else //sort using windows sort
+                       string command = "sort " + outputFile + " /O " + outfile;
+                       system(command.c_str());
+               #endif
+               
+
+               //output to new file distance for each row and save positions in file where new row begins
+               ifstream in;
+               openInputFile(outfile, in);
+               
+               ofstream out;
+               openOutputFile(outputFile, out);
+               
+               out.setf(ios::fixed, ios::floatfield); out.setf(ios::showpoint);
+
+               out << alignDB.getNumSeqs() << endl;
+               
+               //get first currentRow
+               string first, currentRow, second;
+               float dist;
+               int i, j;
+               i = 0; j = 0;
+               map<string, float> rowDists; //take advantage of the fact that maps are already sorted by key 
+               map<string, float>::iterator it;
+               
+               in >> first;
+               currentRow = first;
+               
+               rowDists[first] = 0.00; //distance to yourself is 0.0
+               
+               in.seekg(0);
+               //openInputFile(outfile, in);
+               
+               while(!in.eof()) {
+                       if (m->control_pressed) { in.close(); remove(outfile.c_str()); out.close(); return 0; }
+                       
+                       in >> first >> second >> dist; gobble(in);
+                               
+                       if (first != currentRow) {
+                               //print out last row
+                               out << currentRow << '\t'; //print name
+
+                               //print dists
+                               for (it = rowDists.begin(); it != rowDists.end(); it++) {
+                                       if (j >= i) { break; }
+                                       out << it->second << '\t';
+                                       j++;
+                               }
+                               out << endl;
+                               
+                               //start new row
+                               currentRow = first;
+                               rowDists.clear();
+                               rowDists[first] = 0.00;
+                               rowDists[second] = dist;
+                               j = 0;
+                               i++;
+                       }else{
+                               rowDists[second] = dist;
+                       }
+               }
+               //print out last row
+               out << currentRow << '\t'; //print name
+                               
+               //print dists
+               for (it = rowDists.begin(); it != rowDists.end(); it++) {
+                       out << it->second << '\t';
+               }
+               out << endl;
+               
+               in.close();
+               out.close();
+               
+               remove(outfile.c_str());
+               
+               return 1;
+               
+       }
+       catch(exception& e) {
+               m->errorOut(e, "DistanceCommand", "convertToLowerTriangle");
+               exit(1);
+       }
+}
 /**************************************************************************************************
 void DistanceCommand::appendFiles(string temp, string filename) {
        try{