]> git.donarmstrong.com Git - mothur.git/blobdiff - metastatscommand.cpp
added pcr.seqs command. fixed bug in rarefacftion.single that caused parsing error...
[mothur.git] / metastatscommand.cpp
index 2a55ff08178cb0a7da3d99aff4b157017bb3d27e..4744424d426d61243894d4779fa60dcd20781be3 100644 (file)
@@ -8,9 +8,8 @@
  */
 
 #include "metastatscommand.h"
-#include "metastats.h"
 #include "sharedutilities.h"
-#include "mothurmetastats.h"
+
 
 //**********************************************************************************************************************
 vector<string> MetaStatsCommand::setParameters(){      
@@ -167,7 +166,7 @@ MetaStatsCommand::MetaStatsCommand(string option) {
                        else { 
                                pickedGroups = true;
                                m->splitAtDash(groups, Groups);
-                               m->Groups = Groups;
+                               m->setGroups(Groups);
                        }
                        
                        sets = validParameter.validFile(parameters, "sets", false);                     
@@ -178,14 +177,14 @@ MetaStatsCommand::MetaStatsCommand(string option) {
 
                        
                        string temp = validParameter.validFile(parameters, "iters", false);                     if (temp == "not found") { temp = "1000"; }
-                       convert(temp, iters); 
+                       m->mothurConvert(temp, iters); 
                        
                        temp = validParameter.validFile(parameters, "threshold", false);                        if (temp == "not found") { temp = "0.05"; }
-                       convert(temp, threshold); 
+                       m->mothurConvert(temp, threshold); 
                        
                        temp = validParameter.validFile(parameters, "processors", false);       if (temp == "not found"){       temp = m->getProcessors();      }
                        m->setProcessors(temp);
-                       convert(temp, processors);
+                       m->mothurConvert(temp, processors);
                }
 
        }
@@ -215,8 +214,9 @@ int MetaStatsCommand::execute(){
                //setup the pairwise comparions of sets for metastats
                //calculate number of comparisons i.e. with groups A,B,C = AB, AC, BC = 3;
                //make sure sets are all in designMap
-               SharedUtil* util = new SharedUtil();  
-               util->setGroups(Sets, designMap->namesOfGroups);  
+               SharedUtil* util = new SharedUtil(); 
+               vector<string> dGroups = designMap->getNamesOfGroups();
+               util->setGroups(Sets, dGroups);  
                delete util;
                
                int numGroups = Sets.size();
@@ -231,26 +231,24 @@ int MetaStatsCommand::execute(){
                //only 1 combo
                if (numGroups == 2) { processors = 1; }
                else if (numGroups < 2) { m->mothurOut("Not enough sets, I need at least 2 valid sets. Unable to complete command."); m->mothurOutEndLine(); m->control_pressed = true; }
-               
-               #if defined (__APPLE__) || (__MACH__) || (linux) || (__linux)
-                       if(processors != 1){
-                               int numPairs = namesOfGroupCombos.size();
-                               int numPairsPerProcessor = numPairs / processors;
+
+        if(processors != 1){
+            int numPairs = namesOfGroupCombos.size();
+            int numPairsPerProcessor = numPairs / processors;
                        
-                               for (int i = 0; i < processors; i++) {
-                                       int startPos = i * numPairsPerProcessor;
-                                       if(i == processors - 1){
-                                               numPairsPerProcessor = numPairs - i * numPairsPerProcessor;
-                                       }
-                                       lines.push_back(linePair(startPos, numPairsPerProcessor));
-                               }
-                       }
-               #endif
+            for (int i = 0; i < processors; i++) {
+                int startPos = i * numPairsPerProcessor;
+                if(i == processors - 1){
+                    numPairsPerProcessor = numPairs - i * numPairsPerProcessor;
+                }
+                lines.push_back(linePair(startPos, numPairsPerProcessor));
+            }
+        }
                
                //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) {  outputTypes.clear(); for (int i = 0; i < lookup.size(); i++) {  delete lookup[i];  } m->Groups.clear(); delete input; delete designMap;  for (int i = 0; i < outputNames.size(); i++) {      m->mothurRemove(outputNames[i]); } return 0; }
+                       if (m->control_pressed) {  outputTypes.clear(); for (int i = 0; i < lookup.size(); i++) {  delete lookup[i];  } m->clearGroups(); delete input; delete designMap;  for (int i = 0; i < outputNames.size(); i++) {       m->mothurRemove(outputNames[i]); } return 0; }
        
                        if(allLines == 1 || labels.count(lookup[0]->getLabel()) == 1){                  
 
@@ -281,13 +279,13 @@ int MetaStatsCommand::execute(){
                        //prevent memory leak
                        for (int i = 0; i < lookup.size(); i++) {  delete lookup[i]; lookup[i] = NULL; }
                        
-                       if (m->control_pressed) {  outputTypes.clear(); m->Groups.clear(); delete input;  delete designMap;  for (int i = 0; i < outputNames.size(); i++) {     m->mothurRemove(outputNames[i]); } return 0; }
+                       if (m->control_pressed) {  outputTypes.clear(); m->clearGroups(); delete input;  delete designMap;  for (int i = 0; i < outputNames.size(); i++) {      m->mothurRemove(outputNames[i]); } return 0; }
 
                        //get next line to process
                        lookup = input->getSharedRAbundVectors();                               
                }
                
-               if (m->control_pressed) { outputTypes.clear(); m->Groups.clear(); delete input; delete designMap;  for (int i = 0; i < outputNames.size(); i++) {       m->mothurRemove(outputNames[i]); }  return 0; }
+               if (m->control_pressed) { outputTypes.clear(); m->clearGroups(); delete input; delete designMap;  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;
@@ -315,7 +313,7 @@ int MetaStatsCommand::execute(){
                }
        
                //reset groups parameter
-               m->Groups.clear();  
+               m->clearGroups();  
                delete input; 
                delete designMap;
                
@@ -338,13 +336,13 @@ int MetaStatsCommand::execute(){
 int MetaStatsCommand::process(vector<SharedRAbundVector*>& thisLookUp){
        try {
                
-               #if defined (__APPLE__) || (__MACH__) || (linux) || (__linux)
+               
                                if(processors == 1){
                                        driver(0, namesOfGroupCombos.size(), thisLookUp);
                                }else{
                                        int process = 1;
                                        vector<int> processIDS;
-               
+               #if defined (__APPLE__) || (__MACH__) || (linux) || (__linux) || (__linux__) || (__unix__) || (__unix)
                                        //loop through and create all the processes you want
                                        while (process != processors) {
                                                int pid = fork();
@@ -370,11 +368,67 @@ int MetaStatsCommand::process(vector<SharedRAbundVector*>& thisLookUp){
                                                int temp = processIDS[i];
                                                wait(&temp);
                                        }
-                               }
-               #else
-                               driver(0, namesOfGroupCombos.size(), thisLookUp);
-               #endif
+        #else
+                    
+                    //////////////////////////////////////////////////////////////////////////////////////////////////////
+                    //Windows version shared memory, so be careful when passing variables through the summarySharedData struct. 
+                    //Above fork() will clone, so memory is separate, but that's not the case with windows, 
+                    //Taking advantage of shared memory to pass results vectors.
+                    //////////////////////////////////////////////////////////////////////////////////////////////////////
+                    
+                    vector<metastatsData*> pDataArray; 
+                    DWORD   dwThreadIdArray[processors-1];
+                    HANDLE  hThreadArray[processors-1]; 
+                    
+                    //Create processor worker threads.
+                    for( int i=1; i<processors; i++ ){
+                        
+                        //make copy of lookup so we don't get access violations
+                        vector<SharedRAbundVector*> newLookup;
+                        vector<string> designMapGroups;
+                        for (int k = 0; k < thisLookUp.size(); k++) {
+                            SharedRAbundVector* temp = new SharedRAbundVector();
+                            temp->setLabel(thisLookUp[k]->getLabel());
+                            temp->setGroup(thisLookUp[k]->getGroup());
+                            newLookup.push_back(temp);
+                            designMapGroups.push_back(designMap->getGroup(thisLookUp[k]->getGroup()));
+                        }
+                        
+                        //for each bin
+                        for (int k = 0; k < thisLookUp[0]->getNumBins(); k++) {
+                            if (m->control_pressed) { for (int j = 0; j < newLookup.size(); j++) {  delete newLookup[j];  } return 0; }
+                            for (int j = 0; j < thisLookUp.size(); j++) { newLookup[j]->push_back(thisLookUp[j]->getAbundance(k), thisLookUp[j]->getGroup()); }
+                        }
+                        
+                        // Allocate memory for thread data.
+                        metastatsData* tempSum = new metastatsData(sharedfile, outputDir, m, lines[i].start, lines[i].num, namesOfGroupCombos, newLookup, designMapGroups, iters, threshold);
+                        pDataArray.push_back(tempSum);
+                        processIDS.push_back(i);
+                        
+                        hThreadArray[i-1] = CreateThread(NULL, 0, MyMetastatsThreadFunction, pDataArray[i-1], 0, &dwThreadIdArray[i-1]);   
+                    }
+                    
+                    //do my part
+                                       driver(lines[0].start, lines[0].num, thisLookUp);
+                    
+                    //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++){
+                        for (int j = 0; j < pDataArray[i]->thisLookUp.size(); j++) {  delete pDataArray[i]->thisLookUp[j];  } 
+                        for (int j = 0; j < pDataArray[i]->outputNames.size(); j++) {  
+                            outputNames.push_back(pDataArray[i]->outputNames[j]);
+                            outputTypes["metastats"].push_back(pDataArray[i]->outputNames[j]);
+                        }
+                                                
+                        CloseHandle(hThreadArray[i]);
+                        delete pDataArray[i];
+                    }
+        #endif
 
+                               }
+               
                return 0;
                
        }
@@ -393,17 +447,17 @@ int MetaStatsCommand::driver(int start, int num, vector<SharedRAbundVector*>& th
                        //get set names
                        string setA = namesOfGroupCombos[c][0]; 
                        string setB = namesOfGroupCombos[c][1];
-                       
+               
                        //get filename
                        string outputFileName = outputDir +  m->getRootName(m->getSimpleName(sharedfile)) + thisLookUp[0]->getLabel() + "." + setA + "-" + setB + ".metastats";
                        outputNames.push_back(outputFileName); outputTypes["metastats"].push_back(outputFileName);
-                       int nameLength = outputFileName.length();
-                       char * output = new char[nameLength];
-                       strcpy(output, outputFileName.c_str());
+                       //int nameLength = outputFileName.length();
+                       //char * output = new char[nameLength];
+                       //strcpy(output, outputFileName.c_str());
        
                        //build matrix from shared rabunds
-                       double** data;
-                       data = new double*[thisLookUp[0]->getNumBins()];
+                       //double** data;
+                       //data = new double*[thisLookUp[0]->getNumBins()];
                        
                        vector< vector<double> > data2; data2.resize(thisLookUp[0]->getNumBins());
                        
@@ -428,31 +482,40 @@ int MetaStatsCommand::driver(int start, int num, vector<SharedRAbundVector*>& th
                                m->mothurOut("Missing shared info for " + setA + " or " + setB + ". Skipping comparison."); m->mothurOutEndLine(); 
                                outputNames.pop_back();
                        }else {
+                
+                ofstream outTemp;
+                string tempOut = outputDir + "data." + setA + "-" + setB + ".matrix";
+                m->openOutputFile(tempOut, outTemp);
+                for (int i = 0; i < subset.size(); i++) { outTemp << '\t' << subset[i]->getGroup(); }
+                outTemp << endl;
+                
+                
                                //fill data
                                for (int j = 0; j < thisLookUp[0]->getNumBins(); j++) {
-                                       data[j] = new double[subset.size()];
+                                       //data[j] = new double[subset.size()];
                                        data2[j].resize(subset.size(), 0.0);
+                    outTemp << "OTU" << (j+1);
                                        for (int i = 0; i < subset.size(); i++) {
-                                               data[j][i] = (subset[i]->getAbundance(j));
                                                data2[j][i] = (subset[i]->getAbundance(j));
+                        outTemp << '\t' << subset[i]->getAbundance(j);
                                        }
+                    outTemp << endl;
                                }
-                               
+                               outTemp.close();
                                m->mothurOut("Comparing " + setA + " and " + setB + "..."); m->mothurOutEndLine(); 
-                               metastat_main(output, thisLookUp[0]->getNumBins(), subset.size(), threshold, iters, data, setACount);
+                               //metastat_main(output, thisLookUp[0]->getNumBins(), subset.size(), threshold, iters, data, setACount);
                                
                                m->mothurOutEndLine();
                                MothurMetastats mothurMeta(threshold, iters);
-                               mothurMeta.runMetastats(outputFileName+".myVersion" , data2, setACount);
+                               mothurMeta.runMetastats(outputFileName , data2, setACount);
                                m->mothurOutEndLine();
-                               
                                m->mothurOutEndLine(); 
                        }
                        
                        //free memory
-                       delete output;
-                       for(int i = 0; i < thisLookUp[0]->getNumBins(); i++)  {  delete[] data[i];  }
-                       delete[] data; 
+                       //delete output;
+                       //for(int i = 0; i < thisLookUp[0]->getNumBins(); i++)  {  delete[] data[i];  }
+                       //delete[] data; 
                }
                
                return 0;