]> git.donarmstrong.com Git - mothur.git/blobdiff - filterseqscommand.cpp
added mpi code to cluster.split command
[mothur.git] / filterseqscommand.cpp
index 28e137376c51261498a4d616ae929ba33b1ff3be..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,23 +259,14 @@ int FilterSeqsCommand::filterSequences() {
                                int outMode=MPI_MODE_CREATE|MPI_MODE_WRONLY; 
                                int inMode=MPI_MODE_RDONLY; 
                                
-                               //char* outFilename = new char[filteredFasta.length()];
-                               //memcpy(outFilename, filteredFasta.c_str(), filteredFasta.length());
-                               
                                char outFilename[1024];
                                strcpy(outFilename, filteredFasta.c_str());
-
-                               //char* inFileName = new char[fastafileNames[s].length()];
-                               //memcpy(inFileName, fastafileNames[s].c_str(), 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);
-                               
-                               //delete inFileName;
-                               //delete outFilename;
 
                                if (m->control_pressed) {  MPI_File_close(&inMPI);  MPI_File_close(&outMPI);  return 0;  }
 
@@ -282,8 +276,10 @@ 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;
@@ -303,10 +299,10 @@ 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;
@@ -328,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]);
@@ -356,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; }
@@ -539,7 +542,7 @@ string FilterSeqsCommand::createFilter() {
                
                F.setLength(alignmentLength);
                
-               if(soft != 0 || isTrue(vertical)){
+               if(trump != '*' || isTrue(vertical) || soft != 0){
                        F.initialize();
                }
                
@@ -577,8 +580,10 @@ 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;
@@ -592,11 +597,10 @@ 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;
@@ -611,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;
@@ -633,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;
@@ -664,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{
@@ -697,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);
@@ -708,10 +721,14 @@ 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 = 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;
@@ -721,8 +738,7 @@ string FilterSeqsCommand::createFilter() {
        
        MPI_Barrier(MPI_COMM_WORLD);
 #endif
-               
-               
+                               
                return filterString;
        }
        catch(exception& e) {