]> git.donarmstrong.com Git - mothur.git/blobdiff - distancecommand.cpp
added groups option to get.oturep command
[mothur.git] / distancecommand.cpp
index 23f2834648bb905ff2139c4bc05de59760e4729b..503f64b21cfd0ae64ec7413c7bfc3dc993e7e2d3 100644 (file)
@@ -162,6 +162,8 @@ int DistanceCommand::execute(){
                
                if (abort == true) { return 0; }
                
+               int startTime = time(NULL);
+               
                int numSeqs = alignDB.getNumSeqs();
                cutoff += 0.005;
                
@@ -194,46 +196,111 @@ int DistanceCommand::execute(){
                start = int (sqrt(float(pid)/float(processors)) * numSeqs);
                end = int (sqrt(float(pid+1)/float(processors)) * numSeqs);
                
-               if (pid == 0) { //you are the root process 
-                       //do your part
-                       string outputMyPart;
-                       driverMPI(start, end, outputMyPart, cutoff);
+               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());
                        
-                       ofstream out;
-                       openOutputFile(outputFile, out);
+                       char filename[1024];
+                       strcpy(filename, outputFile.c_str());
                        
-                       out << outputMyPart;
+                       MPI_File_open(MPI_COMM_WORLD, filename, amode, MPI_INFO_NULL, &outMPI);
+                       //delete filename;
+
+                       if (pid == 0) { //you are the root process 
                        
-                       //get the childrens parts
-                       for(int i = 1; i < processors; i++) { 
-                               int length;
-                               MPI_Recv(&length, 1, MPI_INT, i, tag, MPI_COMM_WORLD, &status); 
+                               //do your part
+                               string outputMyPart;
+                               driverMPI(start, end, outMPI, cutoff);
                                
-                               char buf[length];
+                               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; }
                                        
-                               MPI_Recv(buf, length, MPI_CHAR, i, tag, MPI_COMM_WORLD, &status); 
+                                       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);
                                
-                               outputMyPart = buf;
-                               out << outputMyPart;
+                               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);
                        }
                        
-                       out.close();
+                       MPI_File_close(&outMPI);
                        
-               }else { //you are a child process
-                       //do your part
-                       string outputMyPart;
-                       driverMPI(start, end, outputMyPart, cutoff);
-               
-                       //send results to parent
-                       int length = outputMyPart.length();
-                       char buf[length];
-                       strcpy(buf, outputMyPart.c_str()); 
+               }else { //lower triangle format
+                       if (pid == 0) { //you are the root process 
                        
-                       MPI_Send( &length, 1, MPI_INT, 0, tag, MPI_COMM_WORLD);
-                       MPI_Send(buf, length, MPI_CHAR, 0, tag, MPI_COMM_WORLD);
-               }
-               
+                               //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);
+                       }
+               }
+               MPI_Barrier(MPI_COMM_WORLD); //make everyone wait - just in case
 #else          
                                
        #if defined (__APPLE__) || (__MACH__) || (linux) || (__linux)
@@ -276,6 +343,14 @@ int DistanceCommand::execute(){
                
                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
@@ -288,7 +363,7 @@ int DistanceCommand::execute(){
                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;
                
        }
@@ -389,26 +464,18 @@ 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, string& outputString, float cutoff){
+int DistanceCommand::driverMPI(int startLine, int endLine, MPI_File& outMPI, float cutoff){
        try {
-
+               MPI_Status status;
                int startTime = time(NULL);
                
-               outputString = "";
-                               
-               if((output == "lt") && startLine == 0){ outputString += (toString(alignDB.getNumSeqs()) + '\n');        }
+               string outputString = "";
                
                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;  }
@@ -419,32 +486,115 @@ int DistanceCommand::driverMPI(int startLine, int endLine, string& outputString,
                                if(dist <= cutoff){
                                        if (output == "column") { outputString += (alignDB.get(i).getName() + ' ' + alignDB.get(j).getName() + ' ' + toString(dist) + '\n'); }
                                }
-                               if (output == "lt") {   outputString += (toString(dist) + '\t'); }
                                
                                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 (output == "lt") { outputString += '\n'; }
-                       
                        if(i % 100 == 0){
-                               m->mothurOut(toString(i) + "\t" + toString(time(NULL) - startTime)); m->mothurOutEndLine();
+                               //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();
                
+               //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", "driver");
+               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
 /**************************************************************************************************/
 int DistanceCommand::convertMatrix(string outputFile) {
        try{
@@ -533,6 +683,100 @@ int DistanceCommand::convertMatrix(string outputFile) {
                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{