]> git.donarmstrong.com Git - mothur.git/blobdiff - filterseqscommand.cpp
added mpi code to cluster.split command
[mothur.git] / filterseqscommand.cpp
index b992750f1d4fafa1374c9a8a804758cd71ba9cfc..33790fe97632525cbb6c155ab953da531f770078 100644 (file)
@@ -102,6 +102,9 @@ FilterSeqsCommand::FilterSeqsCommand(string option)  {
                        // ...at some point should added some additional type checking...
                        
                        string temp;
+                       hard = validParameter.validFile(parameters, "hard", true);                              if (hard == "not found") { hard = ""; }
+                       else if (hard == "not open") { hard = ""; abort = true; }       
+
                        temp = validParameter.validFile(parameters, "trump", false);                    if (temp == "not found") { temp = "*"; }
                        trump = temp[0];
                        
@@ -111,13 +114,13 @@ FilterSeqsCommand::FilterSeqsCommand(string option)  {
                        temp = validParameter.validFile(parameters, "processors", false);       if (temp == "not found"){       temp = "1";                             }
                        convert(temp, processors); 
                        
-                       hard = validParameter.validFile(parameters, "hard", true);                              if (hard == "not found") { hard = ""; }
-                       else if (hard == "not open") { abort = true; }  
-                       
-                       vertical = validParameter.validFile(parameters, "vertical", false);             if (vertical == "not found") { vertical = "T"; }
+                       vertical = validParameter.validFile(parameters, "vertical", false);             
+                       if (vertical == "not found") { 
+                               if ((hard == "") && (trump == '*') && (soft == 0)) { vertical = "T"; } //you have not given a hard file or set the trump char.
+                               else { vertical = "F";  }
+                       }
                        
                        numSeqs = 0;
-                       
                }
                
        }
@@ -138,11 +141,11 @@ void FilterSeqsCommand::help(){
                m->mothurOut("For example: fasta=abrecovery.fasta-amazon.fasta \n");
                m->mothurOut("The trump parameter .... The default is ...\n");
                m->mothurOut("The soft parameter .... The default is ....\n");
-               m->mothurOut("The hard parameter .... The default is ....\n");
-               m->mothurOut("The vertical parameter .... The default is T.\n");
+               m->mothurOut("The hard parameter allows you to enter a file containing the filter you want to use.\n");
+               m->mothurOut("The vertical parameter removes columns where all sequences contain a gap character. The default is T.\n");
                m->mothurOut("The filter.seqs command should be in the following format: \n");
-               m->mothurOut("filter.seqs(fasta=yourFastaFile, trump=yourTrump, soft=yourSoft, hard=yourHard, vertical=yourVertical) \n");
-               m->mothurOut("Example filter.seqs(fasta=abrecovery.fasta, trump=..., soft=..., hard=..., vertical=T).\n");
+               m->mothurOut("filter.seqs(fasta=yourFastaFile, trump=yourTrump) \n");
+               m->mothurOut("Example filter.seqs(fasta=abrecovery.fasta, trump=.).\n");
                m->mothurOut("Note: No spaces between parameter labels (i.e. fasta), '=' and parameters (i.e.yourFasta).\n\n");
                
        }
@@ -256,15 +259,15 @@ int FilterSeqsCommand::filterSequences() {
                                int outMode=MPI_MODE_CREATE|MPI_MODE_WRONLY; 
                                int inMode=MPI_MODE_RDONLY; 
                                
-                               char outFilename[filteredFasta.length()];
+                               char outFilename[1024];
                                strcpy(outFilename, filteredFasta.c_str());
-                               
-                               char inFileName[fastafileNames[s].length()];
+                       
+                               char inFileName[1024];
                                strcpy(inFileName, fastafileNames[s].c_str());
                                
                                MPI_File_open(MPI_COMM_WORLD, inFileName, inMode, MPI_INFO_NULL, &inMPI);  //comm, filename, mode, info, filepointer
                                MPI_File_open(MPI_COMM_WORLD, outFilename, outMode, MPI_INFO_NULL, &outMPI);
-                               
+
                                if (m->control_pressed) {  MPI_File_close(&inMPI);  MPI_File_close(&outMPI);  return 0;  }
 
                                if (pid == 0) { //you are the root process 
@@ -273,13 +276,16 @@ int FilterSeqsCommand::filterSequences() {
                                        numSeqs += num;
                                        
                                        //send file positions to all processes
-                                       MPI_Bcast(&num, 1, MPI_INT, 0, MPI_COMM_WORLD);  //send numSeqs
-                                       MPI_Bcast(&MPIPos[0], (num+1), MPI_LONG, 0, MPI_COMM_WORLD); //send file pos    
+                                       for(int i = 1; i < processors; i++) { 
+                                               MPI_Send(&num, 1, MPI_INT, i, tag, MPI_COMM_WORLD);
+                                               MPI_Send(&MPIPos[0], (num+1), MPI_LONG, i, tag, MPI_COMM_WORLD);
+                                       }
                                        
                                        //figure out how many sequences you have to do
                                        numSeqsPerProcessor = num / processors;
-                                       if(pid == (processors - 1)){    numSeqsPerProcessor = num - pid * numSeqsPerProcessor;  }
                                        int startIndex =  pid * numSeqsPerProcessor;
+                                       if(pid == (processors - 1)){    numSeqsPerProcessor = num - pid * numSeqsPerProcessor;  }
+                                       
                                
                                        //do your part
                                        driverMPIRun(startIndex, numSeqsPerProcessor, inMPI, outMPI, MPIPos);
@@ -293,15 +299,16 @@ int FilterSeqsCommand::filterSequences() {
                                        }
                                        
                                }else { //you are a child process
-                                       MPI_Bcast(&num, 1, MPI_INT, 0, MPI_COMM_WORLD); //get numSeqs
-                                       numSeqs += num;
+                                       MPI_Recv(&num, 1, MPI_INT, 0, tag, MPI_COMM_WORLD, &status);
                                        MPIPos.resize(num+1);
-                                       MPI_Bcast(&MPIPos[0], (num+1), MPI_LONG, 0, MPI_COMM_WORLD); //get file positions
+                                       numSeqs += num;
+                                       MPI_Recv(&MPIPos[0], (num+1), MPI_LONG, 0, tag, MPI_COMM_WORLD, &status);
                                        
                                        //figure out how many sequences you have to align
                                        numSeqsPerProcessor = num / processors;
-                                       if(pid == (processors - 1)){    numSeqsPerProcessor = num - pid * numSeqsPerProcessor;  }
                                        int startIndex =  pid * numSeqsPerProcessor;
+                                       if(pid == (processors - 1)){    numSeqsPerProcessor = num - pid * numSeqsPerProcessor;  }
+                                       
                                        
                                        //align your part
                                        driverMPIRun(startIndex, numSeqsPerProcessor, inMPI, outMPI, MPIPos);           
@@ -317,17 +324,21 @@ int FilterSeqsCommand::filterSequences() {
                                
                                MPI_File_close(&outMPI);
                                MPI_File_close(&inMPI);
+                               MPI_Barrier(MPI_COMM_WORLD); //make everyone wait - just in case
                                
 #else
                #if defined (__APPLE__) || (__MACH__) || (linux) || (__linux)
                                if(processors == 1){
                                        ifstream inFASTA;
+                                       int numFastaSeqs;
                                        openInputFile(fastafileNames[s], inFASTA);
-                                       int numFastaSeqs=count(istreambuf_iterator<char>(inFASTA),istreambuf_iterator<char>(), '>');
+                                       getNumSeqs(inFASTA, numFastaSeqs);
                                        inFASTA.close();
                                        
                                        lines.push_back(new linePair(0, numFastaSeqs));
                                        
+                                       numSeqs += numFastaSeqs;
+                                       
                                        driverRunFilter(filter, filteredFasta, fastafileNames[s], lines[0]);
                                }else{
                                        setLines(fastafileNames[s]);
@@ -345,12 +356,15 @@ int FilterSeqsCommand::filterSequences() {
                                if (m->control_pressed) {  return 1; }
                #else
                                ifstream inFASTA;
+                               int numFastaSeqs;
                                openInputFile(fastafileNames[s], inFASTA);
-                               int numFastaSeqs=count(istreambuf_iterator<char>(inFASTA),istreambuf_iterator<char>(), '>');
+                               getNumSeqs(inFASTA, numFastaSeqs);
                                inFASTA.close();
                                        
                                lines.push_back(new linePair(0, numFastaSeqs));
                                
+                               numSeqs += numFastaSeqs;
+                               
                                driverRunFilter(filter, filteredFasta, fastafileNames[s], lines[0]);
 
                                if (m->control_pressed) {  return 1; }
@@ -380,12 +394,13 @@ int FilterSeqsCommand::driverMPIRun(int start, int num, MPI_File& inMPI, MPI_Fil
                
                        //read next sequence
                        int length = MPIPos[start+i+1] - MPIPos[start+i];
-                       char buf4[length];
+                       char* buf4 = new char[length];
                        MPI_File_read_at(inMPI, MPIPos[start+i], buf4, length, MPI_CHAR, &status);
                        
                        string tempBuf = buf4;
                        if (tempBuf.length() > length) { tempBuf = tempBuf.substr(0, length);  }
                        istringstream iss (tempBuf,istringstream::in);
+                       delete buf4;
        
                        Sequence seq(iss);  gobble(iss);
                        
@@ -405,11 +420,12 @@ int FilterSeqsCommand::driverMPIRun(int start, int num, MPI_File& inMPI, MPI_Fil
                                if(count % 10 == 0){ //output to file 
                                        //send results to parent
                                        int length = outputString.length();
-                                       char buf[length];
-                                       strcpy(buf, outputString.c_str()); 
+                                       char* buf = new char[length];
+                                       memcpy(buf, outputString.c_str(), length);
                                
                                        MPI_File_write_shared(outMPI, buf, length, MPI_CHAR, &status);
                                        outputString = "";
+                                       delete buf;
                                }
 
                        }
@@ -420,11 +436,12 @@ int FilterSeqsCommand::driverMPIRun(int start, int num, MPI_File& inMPI, MPI_Fil
                if(outputString != ""){ //output to file 
                        //send results to parent
                        int length = outputString.length();
-                       char buf[length];
-                       strcpy(buf, outputString.c_str()); 
+                       char* buf = new char[length];
+                       memcpy(buf, outputString.c_str(), length);
                        
                        MPI_File_write_shared(outMPI, buf, length, MPI_CHAR, &status);
                        outputString = "";
+                       delete buf;
                }
                
                if((num) % 100 != 0){   cout << (num) << endl;   m->mothurOutJustToLog(toString(num) + "\n");   }
@@ -525,7 +542,7 @@ string FilterSeqsCommand::createFilter() {
                
                F.setLength(alignmentLength);
                
-               if(soft != 0 || isTrue(vertical)){
+               if(trump != '*' || isTrue(vertical) || soft != 0){
                        F.initialize();
                }
                
@@ -548,8 +565,11 @@ string FilterSeqsCommand::createFilter() {
                                MPI_Comm_size(MPI_COMM_WORLD, &processors);
                                MPI_Comm_rank(MPI_COMM_WORLD, &pid); 
                                                        
-                               char* tempFileName = new char(fastafileNames[s].length());
-                               tempFileName = &(fastafileNames[s][0]);
+                               //char* tempFileName = new char(fastafileNames[s].length());
+                               //tempFileName = &(fastafileNames[s][0]);
+                               
+                               char tempFileName[1024];
+                               strcpy(tempFileName, fastafileNames[s].c_str());
                
                                MPI_File_open(MPI_COMM_WORLD, tempFileName, MPI_MODE_RDONLY, MPI_INFO_NULL, &inMPI);  //comm, filename, mode, info, filepointer
                                
@@ -560,13 +580,16 @@ string FilterSeqsCommand::createFilter() {
                                                numSeqs += num;
                                                
                                                //send file positions to all processes
-                                               MPI_Bcast(&num, 1, MPI_INT, 0, MPI_COMM_WORLD);  //send numSeqs
-                                               MPI_Bcast(&MPIPos[0], (num+1), MPI_LONG, 0, MPI_COMM_WORLD); //send file pos    
+                                               for(int i = 1; i < processors; i++) { 
+                                                       MPI_Send(&num, 1, MPI_INT, i, tag, MPI_COMM_WORLD);
+                                                       MPI_Send(&MPIPos[0], (num+1), MPI_LONG, i, tag, MPI_COMM_WORLD);
+                                               }
                                                                
                                                //figure out how many sequences you have to do
                                                numSeqsPerProcessor = num / processors;
-                                               if(pid == (processors - 1)){    numSeqsPerProcessor = num - pid * numSeqsPerProcessor;  }
                                                int startIndex =  pid * numSeqsPerProcessor;
+                                               if(pid == (processors - 1)){    numSeqsPerProcessor = num - pid * numSeqsPerProcessor;  }
+                                               
                                
                                                //do your part
                                                MPICreateFilter(startIndex, numSeqsPerProcessor, F, inMPI, MPIPos);
@@ -574,16 +597,16 @@ string FilterSeqsCommand::createFilter() {
                                                if (m->control_pressed) {  MPI_File_close(&inMPI);  return 0;  }
                                                                                                
                                }else { //i am the child process
-                       
-                                       MPI_Bcast(&num, 1, MPI_INT, 0, MPI_COMM_WORLD); //get numSeqs
+                                       MPI_Recv(&num, 1, MPI_INT, 0, tag, MPI_COMM_WORLD, &status);
                                        MPIPos.resize(num+1);
                                        numSeqs += num;
-                                       MPI_Bcast(&MPIPos[0], (num+1), MPI_LONG, 0, MPI_COMM_WORLD); //get file positions
+                                       MPI_Recv(&MPIPos[0], (num+1), MPI_LONG, 0, tag, MPI_COMM_WORLD, &status);
                                        
                                        //figure out how many sequences you have to align
                                        numSeqsPerProcessor = num / processors;
-                                       if(pid == (processors - 1)){    numSeqsPerProcessor = num - pid * numSeqsPerProcessor;  }
                                        int startIndex =  pid * numSeqsPerProcessor;
+                                       if(pid == (processors - 1)){    numSeqsPerProcessor = num - pid * numSeqsPerProcessor;  }
+                                       
                                        
                                        //do your part
                                        MPICreateFilter(startIndex, numSeqsPerProcessor, F, inMPI,  MPIPos);
@@ -592,13 +615,15 @@ string FilterSeqsCommand::createFilter() {
                                }
                                
                                MPI_File_close(&inMPI);
+                               MPI_Barrier(MPI_COMM_WORLD); //make everyone wait - just in case
                                
 #else
                #if defined (__APPLE__) || (__MACH__) || (linux) || (__linux)
                                if(processors == 1){
                                        ifstream inFASTA;
+                                       int numFastaSeqs;
                                        openInputFile(fastafileNames[s], inFASTA);
-                                       int numFastaSeqs=count(istreambuf_iterator<char>(inFASTA),istreambuf_iterator<char>(), '>');
+                                       getNumSeqs(inFASTA, numFastaSeqs);
                                        inFASTA.close();
                                        
                                        numSeqs += numFastaSeqs;
@@ -614,8 +639,9 @@ string FilterSeqsCommand::createFilter() {
                                if (m->control_pressed) {  return filterString; }
                #else
                                ifstream inFASTA;
+                               int numFastaSeqs;
                                openInputFile(fastafileNames[s], inFASTA);
-                               int numFastaSeqs=count(istreambuf_iterator<char>(inFASTA),istreambuf_iterator<char>(), '>');
+                               getNumSeqs(inFASTA, numFastaSeqs);
                                inFASTA.close();
                                
                                numSeqs += numFastaSeqs;
@@ -645,20 +671,24 @@ string FilterSeqsCommand::createFilter() {
                                vector<int> temp; temp.resize(alignmentLength+1);
                                                                
                                //get the frequencies from the child processes
-                               for(int i = 0; i < ((processors-1)*5); i++) { 
-                                       MPI_Recv(&temp[0], (alignmentLength+1), MPI_INT, MPI_ANY_SOURCE, 2001, MPI_COMM_WORLD, &status); 
-                                       int receiveTag = temp[temp.size()-1];  //child process added a int to the end to indicate what letter count this is for
-                                       
-                                       if (receiveTag == Atag) { //you are recieveing the A frequencies
-                                               for (int k = 0; k < alignmentLength; k++) {             F.a[k] += temp[k];      }
-                                       }else if (receiveTag == Ttag) { //you are recieveing the T frequencies
-                                               for (int k = 0; k < alignmentLength; k++) {             F.t[k] += temp[k];      }
-                                       }else if (receiveTag == Ctag) { //you are recieveing the C frequencies
-                                               for (int k = 0; k < alignmentLength; k++) {             F.c[k] += temp[k];      }
-                                       }else if (receiveTag == Gtag) { //you are recieveing the G frequencies
-                                               for (int k = 0; k < alignmentLength; k++) {             F.g[k] += temp[k];      }
-                                       }else if (receiveTag == Gaptag) { //you are recieveing the gap frequencies
-                                               for (int k = 0; k < alignmentLength; k++) {             F.gap[k] += temp[k];    }
+                               for(int i = 1; i < processors; i++) { 
+                               
+                                       for (int j = 0; j < 5; j++) {
+                                       
+                                               MPI_Recv(&temp[0], (alignmentLength+1), MPI_INT, i, 2001, MPI_COMM_WORLD, &status); 
+                                               int receiveTag = temp[temp.size()-1];  //child process added a int to the end to indicate what letter count this is for
+                                               
+                                               if (receiveTag == Atag) { //you are recieveing the A frequencies
+                                                       for (int k = 0; k < alignmentLength; k++) {             F.a[k] += temp[k];      }
+                                               }else if (receiveTag == Ttag) { //you are recieveing the T frequencies
+                                                       for (int k = 0; k < alignmentLength; k++) {             F.t[k] += temp[k];      }
+                                               }else if (receiveTag == Ctag) { //you are recieveing the C frequencies
+                                                       for (int k = 0; k < alignmentLength; k++) {             F.c[k] += temp[k];      }
+                                               }else if (receiveTag == Gtag) { //you are recieveing the G frequencies
+                                                       for (int k = 0; k < alignmentLength; k++) {             F.g[k] += temp[k];      }
+                                               }else if (receiveTag == Gaptag) { //you are recieveing the gap frequencies
+                                                       for (int k = 0; k < alignmentLength; k++) {             F.gap[k] += temp[k];    }
+                                               }
                                        }
                                } 
                        }else{
@@ -678,6 +708,8 @@ string FilterSeqsCommand::createFilter() {
                        
                }
                
+               MPI_Barrier(MPI_COMM_WORLD); //make everyone wait - just in case
+               
                if (pid == 0) { //only one process should output the filter
 #endif
                F.setNumSeqs(numSeqs);
@@ -689,20 +721,24 @@ string FilterSeqsCommand::createFilter() {
                
 #ifdef USE_MPI
                //send filter string to kids
-               MPI_Bcast(&filterString[0], alignmentLength, MPI_CHAR, 0, MPI_COMM_WORLD); 
+               //for(int i = 1; i < processors; i++) { 
+               //      MPI_Send(&filterString[0], alignmentLength, MPI_CHAR, i, 2001, MPI_COMM_WORLD);
+               //}
+               MPI_Bcast(&filterString[0], alignmentLength, MPI_CHAR, 0, MPI_COMM_WORLD);
        }else{
                //recieve filterString
-               char tempBuf[alignmentLength];
+               char* tempBuf = new char[alignmentLength];
+               //MPI_Recv(&tempBuf[0], alignmentLength, MPI_CHAR, 0, 2001, MPI_COMM_WORLD, &status);
                MPI_Bcast(tempBuf, alignmentLength, MPI_CHAR, 0, MPI_COMM_WORLD);
                
                filterString = tempBuf;
                if (filterString.length() > alignmentLength) { filterString = filterString.substr(0, alignmentLength);  }
+               delete tempBuf; 
        }
        
        MPI_Barrier(MPI_COMM_WORLD);
 #endif
-               
-               
+                               
                return filterString;
        }
        catch(exception& e) {
@@ -764,13 +800,14 @@ int FilterSeqsCommand::MPICreateFilter(int start, int num, Filters& F, MPI_File&
                        //read next sequence
                        int length = MPIPos[start+i+1] - MPIPos[start+i];
        
-                       char buf4[length];
+                       char* buf4 = new char[length];
                        MPI_File_read_at(inMPI, MPIPos[start+i], buf4, length, MPI_CHAR, &status);
                        
                        string tempBuf = buf4;
                        if (tempBuf.length() > length) { tempBuf = tempBuf.substr(0, length);  }
                        istringstream iss (tempBuf,istringstream::in);
-       
+                       delete buf4;
+
                        Sequence seq(iss);  
 
                        if (seq.getAligned().length() != alignmentLength) {  cout << "Alignment length is " << alignmentLength << " and sequence " << seq.getName() << " has length " << seq.getAligned().length() << ", please correct." << endl; exit(1);  }