]> git.donarmstrong.com Git - mothur.git/blobdiff - countseqscommand.cpp
added shared file to make.table
[mothur.git] / countseqscommand.cpp
index dfa012eeff0f0f49c778c00b2098348a489f1855..301aff7007c309d18952f04493188b9330b56f55 100644 (file)
@@ -8,15 +8,17 @@
  */
 
 #include "countseqscommand.h"
-#include "groupmap.h"
 #include "sharedutilities.h"
 #include "counttable.h"
+#include "inputdata.h"
 
 //**********************************************************************************************************************
 vector<string> CountSeqsCommand::setParameters(){      
        try {
-               CommandParameter pname("name", "InputTypes", "", "", "none", "none", "none","count",false,true,true); parameters.push_back(pname);
-               CommandParameter pgroup("group", "InputTypes", "", "", "none", "none", "none","",false,false,true); parameters.push_back(pgroup);
+        CommandParameter pshared("shared", "InputTypes", "", "", "NameSHared-sharedGroup", "NameSHared", "none","count",false,false,true); parameters.push_back(pshared);
+               CommandParameter pname("name", "InputTypes", "", "", "NameSHared", "NameSHared", "none","count",false,false,true); parameters.push_back(pname);
+               CommandParameter pgroup("group", "InputTypes", "", "", "sharedGroup", "none", "none","",false,false,true); parameters.push_back(pgroup);
+        CommandParameter pprocessors("processors", "Number", "", "1", "", "", "","",false,false,true); parameters.push_back(pprocessors);
         CommandParameter plarge("large", "Boolean", "", "F", "", "", "","",false,false); parameters.push_back(plarge);
                CommandParameter pgroups("groups", "String", "", "", "", "", "","",false,false); parameters.push_back(pgroups);
                CommandParameter pinputdir("inputdir", "String", "", "", "", "", "","",false,false); parameters.push_back(pinputdir);
@@ -35,10 +37,11 @@ vector<string> CountSeqsCommand::setParameters(){
 string CountSeqsCommand::getHelpString(){      
        try {
                string helpString = "";
-               helpString += "The count.seqs aka. make.table command reads a name file and outputs a .count_table file.  You may also provide a group file to get the counts broken down by group.\n";
+               helpString += "The count.seqs aka. make.table command reads a name or shared file and outputs a .count_table file.  You may also provide a group with the names file to get the counts broken down by group.\n";
                helpString += "The groups parameter allows you to indicate which groups you want to include in the counts, by default all groups in your groupfile are used.\n";
         helpString += "The large parameter indicates the name and group files are too large to fit in RAM.\n";
                helpString += "When you use the groups parameter and a sequence does not represent any sequences from the groups you specify it is not included in the .count.summary file.\n";
+        helpString += "The processors parameter allows you to specify the number of processors to use. The default is 1.\n";
                helpString += "The count.seqs command should be in the following format: count.seqs(name=yourNameFile).\n";
                helpString += "Example count.seqs(name=amazon.names) or make.table(name=amazon.names).\n";
                helpString += "Note: No spaces between parameter labels (i.e. name), '=' and parameters (i.e.yourNameFile).\n";
@@ -53,8 +56,7 @@ string CountSeqsCommand::getHelpString(){
 string CountSeqsCommand::getOutputPattern(string type) {
     try {
         string pattern = "";
-        
-        if (type == "count") {  pattern = "[filename],count_table"; }
+        if (type == "count") {  pattern = "[filename],count_table-[filename],[distance],count_table"; }
         else { m->mothurOut("[ERROR]: No definition for type " + type + " output pattern.\n"); m->control_pressed = true;  }
         
         return pattern;
@@ -81,7 +83,8 @@ CountSeqsCommand::CountSeqsCommand(){
 
 CountSeqsCommand::CountSeqsCommand(string option)  {
        try {
-               abort = false; calledHelp = false;   
+               abort = false; calledHelp = false;
+        allLines = 1;
                
                //allow user to run help
                if(option == "help") { help(); abort = true; calledHelp = true; }
@@ -125,31 +128,58 @@ CountSeqsCommand::CountSeqsCommand(string option)  {
                                        //if the user has not given a path then, add inputdir. else leave path alone.
                                        if (path == "") {       parameters["group"] = inputDir + it->second;            }
                                }
+                
+                it = parameters.find("shared");
+                               //user has given a template file
+                               if(it != parameters.end()){
+                                       path = m->hasPath(it->second);
+                                       //if the user has not given a path then, add inputdir. else leave path alone.
+                                       if (path == "") {       parameters["shared"] = inputDir + it->second;           }
+                               }
                        }
                        
                        //check for required parameters
                        namefile = validParameter.validFile(parameters, "name", true);
                        if (namefile == "not open") { namefile = ""; abort = true; }
-                       else if (namefile == "not found"){                                      
-                               namefile = m->getNameFile(); 
-                               if (namefile != "") { m->mothurOut("Using " + namefile + " as input file for the name parameter."); m->mothurOutEndLine(); }
-                               else {  m->mothurOut("You have no current namefile and the name parameter is required."); m->mothurOutEndLine(); abort = true; }
-                       }else { m->setNameFile(namefile); }
-                       
+                       else if (namefile == "not found"){      namefile = ""; }
+            else { m->setNameFile(namefile); }
+            
+            sharedfile = validParameter.validFile(parameters, "shared", true);
+                       if (sharedfile == "not open") { sharedfile = ""; abort = true; }
+                       else if (sharedfile == "not found"){    sharedfile = ""; }
+            else { m->setSharedFile(sharedfile); }
+            
                        groupfile = validParameter.validFile(parameters, "group", true);
                        if (groupfile == "not open") { abort = true; }
                        else if (groupfile == "not found") {  groupfile = "";  }        
                        else { m->setGroupFile(groupfile); }
-                       
+            
+            if ((namefile == "") && (sharedfile == "")) {
+                namefile = m->getNameFile();
+                               if (namefile != "") { m->mothurOut("Using " + namefile + " as input file for the name parameter."); m->mothurOutEndLine(); }
+                               else {
+                    sharedfile = m->getSharedFile();
+                    if (sharedfile != "") { m->mothurOut("Using " + sharedfile + " as input file for the shared parameter."); m->mothurOutEndLine(); }
+                    else {
+                        m->mothurOut("You have no current namefile or sharedfile and the name or shared parameter is required."); m->mothurOutEndLine(); abort = true;
+                    }
+                }
+                       }
+
                        groups = validParameter.validFile(parameters, "groups", false);                 
                        if (groups == "not found") { groups = "all"; }
                        m->splitAtDash(groups, Groups);
+            m->setGroups(Groups);
             
             string temp = validParameter.validFile(parameters, "large", false);                if (temp == "not found") {      temp = "F";     }
                        large = m->isTrue(temp);
+            
+            temp = validParameter.validFile(parameters, "processors", false);  if (temp == "not found"){       temp = m->getProcessors();      }
+                       m->setProcessors(temp);
+                       m->mothurConvert(temp, processors);
                        
                        //if the user changes the output directory command factory will send this info to us in the output parameter 
-                       outputDir = validParameter.validFile(parameters, "outputdir", false);           if (outputDir == "not found"){  outputDir = m->hasPath(namefile);               }
+                       outputDir = validParameter.validFile(parameters, "outputdir", false);           if (outputDir == "not found"){  outputDir = "";         }
 
                }
                
@@ -166,16 +196,108 @@ int CountSeqsCommand::execute(){
                
                if (abort == true) { if (calledHelp) { return 0; }  return 2;   }
                
-        map<string, string> variables; 
-        variables["[filename]"] = outputDir + m->getRootName(m->getSimpleName(namefile));
-               string outputFileName = getOutputFileName("count", variables);
-               
-        int total = 0;
-        if (!large) { total = processSmall(outputFileName); }
-        else { total = processLarge(outputFileName);  }
+        map<string, string> variables;
+
+        if (namefile != "") {
+            int total = 0;
+            int start = time(NULL);
+            if (outputDir == "") { outputDir = m->hasPath(namefile); }
+            variables["[filename]"] = outputDir + m->getRootName(m->getSimpleName(namefile));
+            string outputFileName = getOutputFileName("count", variables);
+            
+            if (!large) { total = processSmall(outputFileName); }
+            else { total = processLarge(outputFileName);  }
+            
+            if (m->control_pressed) { m->mothurRemove(outputFileName); return 0; }
+            
+            m->mothurOut("It took " + toString(time(NULL) - start) + " secs to create a table for " + toString(total) + " sequences.");
+            m->mothurOutEndLine(); m->mothurOutEndLine();
+            
+            m->mothurOutEndLine();
+            m->mothurOut("Total number of sequences: " + toString(total)); m->mothurOutEndLine();
+        }else {
+            if (outputDir == "") { outputDir = m->hasPath(sharedfile); }
+            variables["[filename]"] = outputDir + m->getRootName(m->getSimpleName(sharedfile));
+            
+            InputData input(sharedfile, "sharedfile");
+            vector<SharedRAbundVector*> lookup = input.getSharedRAbundVectors();
+            string lastLabel = lookup[0]->getLabel();
+            
+            //if the users enters label "0.06" and there is no "0.06" in their file use the next lowest label.
+            set<string> processedLabels;
+            set<string> userLabels = labels;
+            
+            //as long as you are not at the end of the file or done wih the lines you want
+            while((lookup[0] != NULL) && ((allLines == 1) || (userLabels.size() != 0))) {
+                
+                if (m->control_pressed) { for (int i = 0; i < lookup.size(); i++) {  delete lookup[i];  } for (int i = 0; i < outputNames.size(); i++) {       m->mothurRemove(outputNames[i]); } return 0; }
+                
+                if(allLines == 1 || labels.count(lookup[0]->getLabel()) == 1){
+                    
+                    m->mothurOut(lookup[0]->getLabel()); m->mothurOutEndLine();
+                    
+                    processShared(lookup, variables);
+                    
+                    processedLabels.insert(lookup[0]->getLabel());
+                    userLabels.erase(lookup[0]->getLabel());
+                }
+                
+                if ((m->anyLabelsToProcess(lookup[0]->getLabel(), userLabels, "") == true) && (processedLabels.count(lastLabel) != 1)) {
+                    string saveLabel = lookup[0]->getLabel();
+                    
+                    for (int i = 0; i < lookup.size(); i++) {  delete lookup[i];  }
+                    lookup = input.getSharedRAbundVectors(lastLabel);
+                    m->mothurOut(lookup[0]->getLabel()); m->mothurOutEndLine();
+                    
+                    processShared(lookup, variables);
+                    
+                    processedLabels.insert(lookup[0]->getLabel());
+                    userLabels.erase(lookup[0]->getLabel());
+                    
+                    //restore real lastlabel to save below
+                    lookup[0]->setLabel(saveLabel);
+                }
+                
+                lastLabel = lookup[0]->getLabel();
+                //prevent memory leak
+                for (int i = 0; i < lookup.size(); i++) {  delete lookup[i]; lookup[i] = NULL; }
+                
+                if (m->control_pressed) { return 0; }
+                
+                //get next line to process
+                lookup = input.getSharedRAbundVectors();
+            }
+            
+            if (m->control_pressed) { for (int i = 0; i < outputNames.size(); i++) {   m->mothurRemove(outputNames[i]); }  return 0; }
+            
+            //output error messages about any remaining user labels
+            set<string>::iterator it;
+            bool needToRun = false;
+            for (it = userLabels.begin(); it != userLabels.end(); it++) {
+                m->mothurOut("Your file does not include the label " + *it);
+                if (processedLabels.count(lastLabel) != 1) {
+                    m->mothurOut(". I will use " + lastLabel + "."); m->mothurOutEndLine();
+                    needToRun = true;
+                }else {
+                    m->mothurOut(". Please refer to " + lastLabel + "."); m->mothurOutEndLine();
+                }
+            }
+            
+            //run last label if you need to
+            if (needToRun == true)  {
+                for (int i = 0; i < lookup.size(); i++) { if (lookup[i] != NULL) { delete lookup[i]; } }
+                lookup = input.getSharedRAbundVectors(lastLabel);
+                
+                m->mothurOut(lookup[0]->getLabel()); m->mothurOutEndLine();
+                
+                processShared(lookup, variables);
+                
+                for (int i = 0; i < lookup.size(); i++) {  delete lookup[i];  }
+            }
+            
+        }
         
-               if (m->control_pressed) { m->mothurRemove(outputFileName); return 0; }
-               
         //set rabund file as new current rabundfile
                itTypes = outputTypes.find("count");
                if (itTypes != outputTypes.end()) {
@@ -183,12 +305,10 @@ int CountSeqsCommand::execute(){
                }
         
         m->mothurOutEndLine();
-               m->mothurOut("Total number of sequences: " + toString(total)); m->mothurOutEndLine();
-               m->mothurOutEndLine();
                m->mothurOut("Output File Names: "); m->mothurOutEndLine();
-               m->mothurOut(outputFileName); m->mothurOutEndLine();    
+               for(int i = 0; i < outputNames.size(); i++) {  m->mothurOut(outputNames[i]); m->mothurOutEndLine();      }
                m->mothurOutEndLine();
-               
+        
                return 0;               
        }
        
@@ -199,6 +319,40 @@ int CountSeqsCommand::execute(){
 }
 //**********************************************************************************************************************
 
+int CountSeqsCommand::processShared(vector<SharedRAbundVector*>& lookup, map<string, string> variables){
+    try {
+        variables["[distance]"] = lookup[0]->getLabel();
+        string outputFileName = getOutputFileName("count", variables);
+        outputNames.push_back(outputFileName); outputTypes["count"].push_back(outputFileName);
+        
+        ofstream out;
+        m->openOutputFile(outputFileName, out);
+        
+        out << "OTU_Label\ttotal\t";
+        for (int i = 0; i < lookup.size(); i++) { out << lookup[i]->getGroup() << '\t'; } out << endl;
+        
+        for (int j = 0; j < lookup[0]->getNumBins(); j++) {
+            if (m->control_pressed) { break; }
+            
+            int total = 0;
+            string output = "";
+            for (int i = 0; i < lookup.size(); i++) {
+                total += lookup[i]->getAbundance(j);
+                output += toString(lookup[i]->getAbundance(j)) + '\t';
+            }
+            out << m->currentSharedBinLabels[j] << '\t' << total << '\t' << output << endl;
+        }
+        out.close();
+        
+        return 0;
+    }
+    catch(exception& e) {
+        m->errorOut(e, "CountSeqsCommand", "processShared");
+        exit(1);
+    }
+}
+//**********************************************************************************************************************
+
 int CountSeqsCommand::processSmall(string outputFileName){
        try {
         ofstream out;
@@ -225,13 +379,172 @@ int CountSeqsCommand::processSmall(string outputFileName){
                        }
                }
                out << endl;
+        out.close();
+        
+        int total = createProcesses(groupMap, outputFileName);
+        
+        if (groupfile != "") { delete groupMap; }
+        
+        return total;
+    }
+       catch(exception& e) {
+               m->errorOut(e, "CountSeqsCommand", "processSmall");
+               exit(1);
+       }
+}
+/**************************************************************************************************/
+int CountSeqsCommand::createProcesses(GroupMap*& groupMap, string outputFileName) {
+       try {
                
-               //open input file
-               ifstream in;
+               vector<int> processIDS;
+               int process = 0;
+        vector<unsigned long long> positions;
+        vector<linePair> lines;
+        int numSeqs = 0;
+        
+#if defined (__APPLE__) || (__MACH__) || (linux) || (__linux) || (__linux__) || (__unix__) || (__unix)
+               positions = m->divideFilePerLine(namefile, processors);
+               for (int i = 0; i < (positions.size()-1); i++) { lines.push_back(linePair(positions[i], positions[(i+1)])); }
+#else
+               if(processors == 1){ lines.push_back(linePair(0, 1000));  }
+        else {
+            int numSeqs = 0;
+            positions = m->setFilePosEachLine(namefile, numSeqs);
+            if (positions.size() < processors) { processors = positions.size(); }
+            
+            //figure out how many sequences you have to process
+            int numSeqsPerProcessor = numSeqs / processors;
+            for (int i = 0; i < processors; i++) {
+                int startIndex =  i * numSeqsPerProcessor;
+                if(i == (processors - 1)){     numSeqsPerProcessor = numSeqs - i * numSeqsPerProcessor;        }
+                lines.push_back(linePair(positions[startIndex], numSeqsPerProcessor));
+            }
+        }
+#endif
+
+                       
+#if defined (__APPLE__) || (__MACH__) || (linux) || (__linux) || (__linux__) || (__unix__) || (__unix)
+               
+               //loop through and create all the processes you want
+               while (process != processors-1) {
+                       int pid = fork();
+                       
+                       if (pid > 0) {
+                               processIDS.push_back(pid);  //create map from line number to pid so you can append files in correct order later
+                               process++;
+                       }else if (pid == 0){
+                string filename = toString(getpid()) + ".temp";
+                               numSeqs = driver(lines[process].start, lines[process].end, filename, groupMap);
+                
+                string tempFile = toString(getpid()) + ".num.temp";
+                ofstream outTemp;
+                m->openOutputFile(tempFile, outTemp);
+                
+                outTemp << numSeqs << endl;
+                outTemp.close();
+                
+                               exit(0);
+                       }else {
+                               m->mothurOut("[ERROR]: unable to spawn the necessary processes."); m->mothurOutEndLine();
+                               for (int i = 0; i < processIDS.size(); i++) { kill (processIDS[i], SIGINT); }
+                               exit(0);
+                       }
+               }
+               
+               string filename = toString(getpid()) + ".temp";
+        numSeqs = driver(lines[processors-1].start, lines[processors-1].end, filename, groupMap);
+        
+               //force parent to wait until all the processes are done
+               for (int i=0;i<processIDS.size();i++) {
+                       int temp = processIDS[i];
+                       wait(&temp);
+               }
+        
+        for (int i = 0; i < processIDS.size(); i++) {
+            string tempFile = toString(processIDS[i]) +  ".num.temp";
+            ifstream intemp;
+            m->openInputFile(tempFile, intemp);
+            
+            int num;
+            intemp >> num; intemp.close();
+            numSeqs += num;
+            m->mothurRemove(tempFile);
+        }
+#else          
+               vector<countData*> pDataArray;
+               DWORD   dwThreadIdArray[processors-1];
+               HANDLE  hThreadArray[processors-1];
+        vector<GroupMap*> copies;
+               
+               //Create processor worker threads.
+               for( int i=0; i<processors-1; i++ ){
+                       string filename = toString(i) + ".temp";
+            
+            GroupMap* copyGroup = new GroupMap();
+            copyGroup->getCopy(groupMap);
+            copies.push_back(copyGroup);
+            vector<string> cGroups = Groups;
+           
+                       countData* temp = new countData(filename, copyGroup, m, lines[i].start, lines[i].end, groupfile, namefile, cGroups);
+                       pDataArray.push_back(temp);
+                       processIDS.push_back(i);
+                       
+                       hThreadArray[i] = CreateThread(NULL, 0, MyCountThreadFunction, pDataArray[i], 0, &dwThreadIdArray[i]);
+               }
+               
+               string filename = toString(processors-1) + ".temp";
+        numSeqs = driver(lines[processors-1].start, lines[processors-1].end, filename, groupMap);
+                       
+               //Wait until all threads have terminated.
+               WaitForMultipleObjects(processors-1, hThreadArray, TRUE, INFINITE);
+               
+               //Close all thread handles and free memory allocations.
+               for(int i=0; i < pDataArray.size(); i++){
+            numSeqs += pDataArray[i]->total;
+            delete copies[i];
+                       CloseHandle(hThreadArray[i]);
+                       delete pDataArray[i];
+               }
+#endif
+               
+               //append output files
+               for(int i=0;i<processIDS.size();i++){
+                       m->appendFiles((toString(processIDS[i]) + ".temp"), outputFileName);
+                       m->mothurRemove((toString(processIDS[i]) + ".temp"));
+               }
+        m->appendFiles(filename, outputFileName);
+        m->mothurRemove(filename);
+
+        
+        //sanity check
+        if (groupfile != "") {
+            if (numSeqs != groupMap->getNumSeqs()) {
+                m->mothurOut("[ERROR]: processes reported processing " + toString(numSeqs) + " sequences, but group file indicates you have " + toString(groupMap->getNumSeqs()) + " sequences.");
+                if (processors == 1) { m->mothurOut(" Could you have a file mismatch?\n"); }
+                else { m->mothurOut(" Either you have a file mismatch or a process failed to complete the task assigned to it.\n"); m->control_pressed = true; }
+            }
+               }
+               return numSeqs;
+       }
+       catch(exception& e) {
+               m->errorOut(e, "CountSeqsCommand", "createProcesses");
+               exit(1);
+       }
+}
+/**************************************************************************************************/
+int CountSeqsCommand::driver(unsigned long long start, unsigned long long end, string outputFileName, GroupMap*& groupMap) {
+       try {
+        
+        ofstream out;
+        m->openOutputFile(outputFileName, out);
+        
+        ifstream in;
                m->openInputFile(namefile, in);
+               in.seekg(start);
         
-               int total = 0;
-               while (!in.eof()) {
+               bool done = false;
+        int total = 0;
+               while (!done) {
                        if (m->control_pressed) { break; }
                        
                        string firstCol, secondCol;
@@ -240,7 +553,7 @@ int CountSeqsCommand::processSmall(string outputFileName){
             m->checkName(firstCol);
             m->checkName(secondCol);
                        //cout << firstCol << '\t' << secondCol << endl;
-           
+            
                        vector<string> names;
                        m->splitAtChar(secondCol, names, ',');
                        
@@ -278,16 +591,23 @@ int CountSeqsCommand::processSmall(string outputFileName){
                        }
                        
                        total += names.size();
+            
+#if defined (__APPLE__) || (__MACH__) || (linux) || (__linux) || (__linux__) || (__unix__) || (__unix)
+            unsigned long long pos = in.tellg();
+            if ((pos == -1) || (pos >= end)) { break; }
+#else
+            if (in.eof()) { break; }
+#endif
+
                }
                in.close();
         out.close();
-               
-               if (groupfile != "") { delete groupMap; }
-
+        
         return total;
+
     }
        catch(exception& e) {
-               m->errorOut(e, "CountSeqsCommand", "processSmall");
+               m->errorOut(e, "CountSeqsCommand", "driver");
                exit(1);
        }
 }
@@ -452,7 +772,7 @@ map<int, string> CountSeqsCommand::processNameFile(string name) {
             }
                }
                in.close();
-        out.close();
+       
                
         if (rest != "") {
             vector<string> pieces = m->splitWhiteSpace(rest);
@@ -475,6 +795,7 @@ map<int, string> CountSeqsCommand::processNameFile(string name) {
             }
 
         }
+        out.close();
         
         return indexToNames;
     }
@@ -528,7 +849,7 @@ map<int, string> CountSeqsCommand::getGroupNames(string filename, set<string>& n
             }
                }
                in.close();
-        out.close();
+        
         
         if (rest != "") {
             vector<string> pieces = m->splitWhiteSpace(rest);
@@ -550,6 +871,7 @@ map<int, string> CountSeqsCommand::getGroupNames(string filename, set<string>& n
                 }
             }
         }
+        out.close();
                
         for (it = groupIndex.begin(); it != groupIndex.end(); it++) {  indexToGroups[it->second] = it->first;  }