changed how we count sequences in a fastafile to allow for '>' in sequence names

[mothur.git] / filterseqscommand.cpp

diff --git a/filterseqscommand.cpp b/filterseqscommand.cpp
index b64e2f530ff91af78cc69dea71b9d2f53f3f981b..a6290a6afe55a28001e5e229e518a7c03c6a1c28 100644 (file)
--- a/filterseqscommand.cpp
+++ b/filterseqscommand.cpp
@@ -282,13 +282,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);
@@ -302,15 +305,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);           
@@ -326,17 +330,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]);
@@ -354,12 +362,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; }
@@ -537,7 +548,7 @@ string FilterSeqsCommand::createFilter() {
                
                F.setLength(alignmentLength);
                
-               if(soft != 0 || isTrue(vertical)){
+               if(trump != '*' || isTrue(vertical) || soft != 0){
                        F.initialize();
                }
                
@@ -575,13 +586,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);
@@ -589,16 +603,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);
@@ -607,13 +621,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;
@@ -629,8 +645,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;
@@ -660,20 +677,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{
@@ -693,6 +714,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);
@@ -704,10 +727,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;
@@ -717,8 +744,7 @@ string FilterSeqsCommand::createFilter() {
        
        MPI_Barrier(MPI_COMM_WORLD);
 #endif
-               
-               
+                               
                return filterString;
        }
        catch(exception& e) {