added root parameter to the unifrac commands so you can choose to include the entire...

[mothur.git] / shhhercommand.cpp

diff --git a/shhhercommand.cpp b/shhhercommand.cpp
index 5bf802ef304dde8c9d628db63a9e0915f23f236e..066a9c51f26b764f14fa3f9490e094748ad3e7b6 100644 (file)
--- a/shhhercommand.cpp
+++ b/shhhercommand.cpp
@@ -32,7 +32,7 @@
 vector<string> ShhherCommand::getValidParameters(){    
        try {
                string Array[] =  {     
-                       "file", "flow", "lookup", "cutoff", "sigma", "outputdir","inputdir", "processors"       
+                       "file", "flow", "lookup", "cutoff", "sigma", "outputdir","inputdir", "processors", "maxiter", "mindelta"        
                };
                
                vector<string> myArray (Array, Array+(sizeof(Array)/sizeof(string)));
@@ -111,7 +111,7 @@ ShhherCommand::ShhherCommand(string option) {
                        
                        //valid paramters for this command
                        string AlignArray[] =  {
-                               "file", "flow", "lookup", "cutoff", "sigma", "outputdir","inputdir", "processors"       
+                               "file", "flow", "lookup", "cutoff", "sigma", "outputdir","inputdir", "processors", "maxiter", "mindelta"        
                        };
                        
                        vector<string> myArray (AlignArray, AlignArray+(sizeof(AlignArray)/sizeof(string)));
@@ -250,9 +250,6 @@ int ShhherCommand::execute(){
                int tag = 1976;
                MPI_Status status; 
 
-               double begClock = clock();
-               unsigned long int begTime = time(NULL);
-               
                if(pid == 0){
 
                        for(int i=1;i<ncpus;i++){
@@ -288,12 +285,15 @@ int ShhherCommand::execute(){
                        }
                        
                        for(int i=0;i<numFiles;i++){
+                               double begClock = clock();
+                               unsigned long int begTime = time(NULL);
+
                                flowFileName = flowFileVector[i];
-                       
+                               
                                m->mothurOut("\n>>>>>\tProcessing " + flowFileName + " (file " + toString(i+1) + " of " + toString(numFiles) + ")\t<<<<<\n");
                                m->mothurOut("Reading flowgrams...\n");
-                               
                                getFlowData();
+
                                m->mothurOut("Identifying unique flowgrams...\n");
                                getUniques();
 
@@ -368,7 +368,7 @@ int ShhherCommand::execute(){
                                                MPI_Send(&nSeqsPerOTU[0], numOTUs, MPI_INT, i, tag, MPI_COMM_WORLD);
                                                MPI_Send(&cumNumSeqs[0], numOTUs, MPI_INT, i, tag, MPI_COMM_WORLD);
                                        }
-                               
+                                       
                                        calcCentroidsDriver(0, numOTUsOnCPU);
                                        
                                        for(int i=1;i<ncpus;i++){
@@ -487,7 +487,6 @@ int ShhherCommand::execute(){
                }
                else{
                        int abort = 1;
-                       bool live = 1;
 
                        MPI_Recv(&abort, 1, MPI_INT, 0, tag, MPI_COMM_WORLD, &status);
                        if(abort){      return 0;       }
@@ -497,6 +496,8 @@ int ShhherCommand::execute(){
 
                        for(int i=0;i<numFiles;i++){
                                //Now into the pyrodist part
+                               bool live = 1;
+
                                char fileName[1024];
                                MPI_Recv(&fileName, 1024, MPI_CHAR, 0, tag, MPI_COMM_WORLD, &status);
                                MPI_Recv(&numSeqs, 1, MPI_INT, 0, tag, MPI_COMM_WORLD, &status);
@@ -552,12 +553,12 @@ int ShhherCommand::execute(){
                                int total;
 
                                while(live){
-
+                                       
                                        MPI_Recv(&total, 1, MPI_INT, 0, tag, MPI_COMM_WORLD, &status);
                                        singleTau.assign(total, 0.0000);
                                        seqNumber.assign(total, 0);
                                        seqIndex.assign(total, 0);
-                                       
+
                                        MPI_Recv(&change[0], numOTUs, MPI_SHORT, 0, tag, MPI_COMM_WORLD, &status);
                                        MPI_Recv(&centroids[0], numOTUs, MPI_INT, 0, tag, MPI_COMM_WORLD, &status);
                                        MPI_Recv(&singleTau[0], total, MPI_DOUBLE, 0, tag, MPI_COMM_WORLD, &status);
@@ -565,13 +566,12 @@ int ShhherCommand::execute(){
                                        MPI_Recv(&seqIndex[0], total, MPI_INT, 0, tag, MPI_COMM_WORLD, &status);
                                        MPI_Recv(&nSeqsPerOTU[0], total, MPI_INT, 0, tag, MPI_COMM_WORLD, &status);
                                        MPI_Recv(&cumNumSeqs[0], numOTUs, MPI_INT, 0, tag, MPI_COMM_WORLD, &status);
-                                       
+
                                        calcCentroidsDriver(startOTU, endOTU);
 
                                        MPI_Send(&centroids[0], numOTUs, MPI_INT, 0, tag, MPI_COMM_WORLD);
                                        MPI_Send(&change[0], numOTUs, MPI_SHORT, 0, tag, MPI_COMM_WORLD);
 
-                                       
                                        MPI_Recv(&centroids[0], numOTUs, MPI_INT, 0, tag, MPI_COMM_WORLD, &status);
                                        MPI_Recv(&weight[0], numOTUs, MPI_DOUBLE, 0, tag, MPI_COMM_WORLD, &status);
                                        MPI_Recv(&change[0], numOTUs, MPI_SHORT, 0, tag, MPI_COMM_WORLD, &status);
@@ -590,8 +590,8 @@ int ShhherCommand::execute(){
                                }
                        }
                }               
-
                MPI_Barrier(MPI_COMM_WORLD);
+
                return 0;
 
        }
@@ -761,6 +761,11 @@ void ShhherCommand::getFlowData(){
                m->openInputFile(flowFileName, flowFile);
                
                string seqName;
+               seqNameVector.clear();
+               lengths.clear();
+               flowDataIntI.clear();
+               nameMap.clear();
+               
                
                int currentNumFlowCells;
                
@@ -943,7 +948,7 @@ void ShhherCommand::getUniques(){
                uniqueFlowDataIntI.resize(numFlowCells * numUniques);
                uniqueLengths.resize(numUniques);       
                
-               flowDataPrI.assign(numSeqs * numFlowCells, 0);
+               flowDataPrI.resize(numSeqs * numFlowCells, 0);
                for(int i=0;i<flowDataPrI.size();i++)   {       flowDataPrI[i] = getProbIntensity(flowDataIntI[i]);             }
        }
        catch(exception& e) {
@@ -1097,7 +1102,7 @@ string ShhherCommand::createDistFile(int processors){
 
                m->mothurOutEndLine();
                
-               m-mothurOut("Total time: " + toString(time(NULL) - begTime) + '\t' + toString((clock() - begClock)/CLOCKS_PER_SEC) + '\n');
+               m->mothurOut("Total time: " + toString(time(NULL) - begTime) + '\t' + toString((clock() - begClock)/CLOCKS_PER_SEC) + '\n');
                
 
                return fDistFileName;
@@ -1203,7 +1208,11 @@ void ShhherCommand::getOTUData(string listFileName){
                otuData.assign(numSeqs, 0);
                cumNumSeqs.assign(numOTUs, 0);
                nSeqsPerOTU.assign(numOTUs, 0);
-               aaP.resize(numOTUs);
+               aaP.clear();aaP.resize(numOTUs);
+               
+               seqNumber.clear();
+               aaI.clear();
+               seqIndex.clear();
                
                string singleOTU = "";
                
@@ -1255,6 +1264,8 @@ void ShhherCommand::getOTUData(string listFileName){
                        for(int j=nSeqsPerOTU[i];j<numSeqs;j++){
                                aaP[i].push_back(0);
                        }
+                       
+                       
                }
                
                for(int i=1;i<numOTUs;i++){
@@ -1264,6 +1275,7 @@ void ShhherCommand::getOTUData(string listFileName){
                seqIndex = seqNumber;
                
                listFile.close();       
+               
        }
        catch(exception& e) {
                m->errorOut(e, "ShhherCommand", "getOTUData");
@@ -1382,6 +1394,7 @@ void ShhherCommand::calcCentroidsDriver(int start, int finish){
        
        try{
                
+       
                for(int i=start;i<finish;i++){
                        
                        double count = 0;
@@ -1393,7 +1406,7 @@ void ShhherCommand::calcCentroidsDriver(int start, int finish){
                        for(int j=0;j<nSeqsPerOTU[i];j++){
                                count += singleTau[seqNumber[cumNumSeqs[i] + j]];
                        }
-                       
+
                        if(nSeqsPerOTU[i] > 0 && count > MIN_COUNT){
                                vector<double> adF(nSeqsPerOTU[i]);
                                vector<int> anL(nSeqsPerOTU[i]);
@@ -1673,9 +1686,9 @@ void ShhherCommand::calcNewDistancesChildMPI(int startSeq, int stopSeq, vector<i
        try{
                vector<double> newTau(numOTUs,0);
                vector<double> norms(numSeqs, 0);
-               otuIndex.resize(0);
-               seqIndex.resize(0);
-               singleTau.resize(0);
+               otuIndex.clear();
+               seqIndex.clear();
+               singleTau.clear();