]> git.donarmstrong.com Git - mothur.git/blobdiff - matrixoutputcommand.cpp
paralellized dist.shared for windows
[mothur.git] / matrixoutputcommand.cpp
index 4ea8f021d3046a3f25be73b201869a73b98898ce..6537963bb58e034ef00da3812cc6af2bb167be71 100644 (file)
@@ -8,47 +8,6 @@
  */
 
 #include "matrixoutputcommand.h"
-#include "sharedsobscollectsummary.h"
-#include "sharedchao1.h"
-#include "sharedace.h"
-#include "sharednseqs.h"
-#include "sharedjabund.h"
-#include "sharedsorabund.h"
-#include "sharedjclass.h"
-#include "sharedsorclass.h"
-#include "sharedjest.h"
-#include "sharedsorest.h"
-#include "sharedthetayc.h"
-#include "sharedthetan.h"
-#include "sharedkstest.h"
-#include "whittaker.h"
-#include "sharedochiai.h"
-#include "sharedanderbergs.h"
-#include "sharedkulczynski.h"
-#include "sharedkulczynskicody.h"
-#include "sharedlennon.h"
-#include "sharedmorisitahorn.h"
-#include "sharedbraycurtis.h"
-#include "sharedjackknife.h"
-#include "whittaker.h"
-#include "odum.h"
-#include "canberra.h"
-#include "structeuclidean.h"
-#include "structchord.h"
-#include "hellinger.h"
-#include "manhattan.h"
-#include "structpearson.h"
-#include "soergel.h"
-#include "spearman.h"
-#include "structkulczynski.h"
-#include "structchi2.h"
-#include "speciesprofile.h"
-#include "hamming.h"
-#include "gower.h"
-#include "memchi2.h"
-#include "memchord.h"
-#include "memeuclidean.h"
-#include "mempearson.h"
 
 //**********************************************************************************************************************
 vector<string> MatrixOutputCommand::setParameters(){   
@@ -58,6 +17,7 @@ vector<string> MatrixOutputCommand::setParameters(){
                CommandParameter pgroups("groups", "String", "", "", "", "", "",false,false); parameters.push_back(pgroups);
                CommandParameter pcalc("calc", "Multiple", "sharedsobs-sharedchao-sharedace-jabund-sorabund-jclass-sorclass-jest-sorest-thetayc-thetan-kstest-sharednseqs-ochiai-anderberg-kulczynski-kulczynskicody-lennon-morisitahorn-braycurtis-whittaker-odum-canberra-structeuclidean-structchord-hellinger-manhattan-structpearson-soergel-spearman-structkulczynski-speciesprofile-hamming-structchi2-gower-memchi2-memchord-memeuclidean-mempearson", "jclass-thetayc", "", "", "",true,false); parameters.push_back(pcalc);
                CommandParameter poutput("output", "Multiple", "lt-square", "lt", "", "", "",false,false); parameters.push_back(poutput);
+               CommandParameter pprocessors("processors", "Number", "", "1", "", "", "",false,false); parameters.push_back(pprocessors);
                CommandParameter pinputdir("inputdir", "String", "", "", "", "", "",false,false); parameters.push_back(pinputdir);
                CommandParameter poutputdir("outputdir", "String", "", "", "", "", "",false,false); parameters.push_back(poutputdir);
                
@@ -75,7 +35,7 @@ string MatrixOutputCommand::getHelpString(){
        try {
                string helpString = "";
                ValidCalculators validCalculator;
-               helpString += "The dist.shared command parameters are shared, groups, calc, output and label.  shared is a required, unless you have a valid current file.\n";
+               helpString += "The dist.shared command parameters are shared, groups, calc, output, processors and label.  shared is a required, unless you have a valid current file.\n";
                helpString += "The groups parameter allows you to specify which of the groups in your groupfile you would like included used.\n";
                helpString += "The group names are separated by dashes. The label parameter allows you to select what distance levels you would like distance matrices created for, and is also separated by dashes.\n";
                helpString += "The dist.shared command should be in the following format: dist.shared(groups=yourGroups, calc=yourCalcs, label=yourLabels).\n";
@@ -156,6 +116,7 @@ MatrixOutputCommand::MatrixOutputCommand(string option)  {
                                if (sharedfile != "") { m->mothurOut("Using " + sharedfile + " as input file for the shared parameter."); m->mothurOutEndLine(); }
                                else {  m->mothurOut("You have no current sharedfile and the shared parameter is required."); m->mothurOutEndLine(); abort = true; }
                        }else if (sharedfile == "not open") { sharedfile = ""; abort = true; }
+                       else { m->setSharedFile(sharedfile); }
                        
                        //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"){  
@@ -179,8 +140,12 @@ MatrixOutputCommand::MatrixOutputCommand(string option)  {
                        if (groups == "not found") { groups = ""; }
                        else { 
                                m->splitAtDash(groups, Groups);
-                               m->Groups = Groups;
+                               m->setGroups(Groups);
                        }
+                       
+                       string temp = validParameter.validFile(parameters, "processors", false);        if (temp == "not found"){       temp = m->getProcessors();      }
+                       m->setProcessors(temp);
+                       m->mothurConvert(temp, processors); 
                                
                        calc = validParameter.validFile(parameters, "calc", false);                     
                        if (calc == "not found") { calc = "jclass-thetayc";  }
@@ -188,6 +153,11 @@ MatrixOutputCommand::MatrixOutputCommand(string option)  {
                                 if (calc == "default")  {  calc = "jclass-thetayc";  }
                        }
                        m->splitAtDash(calc, Estimators);
+                       if (m->inUsersGroups("citation", Estimators)) { 
+                               ValidCalculators validCalc; validCalc.printCitations(Estimators); 
+                               //remove citation from list of calcs
+                               for (int i = 0; i < Estimators.size(); i++) { if (Estimators[i] == "citation") {  Estimators.erase(Estimators.begin()+i); break; } }
+                       }
 
                        if (abort == false) {
                        
@@ -313,13 +283,18 @@ int MatrixOutputCommand::execute(){
                if (lookup.size() < 2) { m->mothurOut("You have not provided enough valid groups.  I cannot run the command."); m->mothurOutEndLine(); delete input; for (int i = 0; i < lookup.size(); i++) {  delete lookup[i];  } return 0;}
                
                numGroups = lookup.size();
+               lines.resize(processors);
+               for (int i = 0; i < processors; i++) {
+                       lines[i].start = int (sqrt(float(i)/float(processors)) * numGroups);
+                       lines[i].end = int (sqrt(float(i+1)/float(processors)) * numGroups);
+               }       
                
-               if (m->control_pressed) { delete input; for (int i = 0; i < lookup.size(); i++) {  delete lookup[i];  } m->Groups.clear(); return 0;  }
+               if (m->control_pressed) { delete input; for (int i = 0; i < lookup.size(); i++) {  delete lookup[i];  } m->clearGroups(); return 0;  }
                                
                //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(); delete input; for (int i = 0; i < lookup.size(); i++) {  delete lookup[i];  } for (int i = 0; i < outputNames.size(); i++) {     remove(outputNames[i].c_str()); } m->Groups.clear(); return 0;  }
+                       if (m->control_pressed) { outputTypes.clear(); delete input; for (int i = 0; i < lookup.size(); i++) {  delete lookup[i];  } for (int i = 0; i < outputNames.size(); i++) {     m->mothurRemove(outputNames[i]); } m->clearGroups(); return 0;  }
                
                        if(allLines == 1 || labels.count(lookup[0]->getLabel()) == 1){                  
                                m->mothurOut(lookup[0]->getLabel()); m->mothurOutEndLine();
@@ -352,7 +327,7 @@ int MatrixOutputCommand::execute(){
                        lookup = input->getSharedRAbundVectors();
                }
                
-               if (m->control_pressed) { outputTypes.clear(); delete input; for (int i = 0; i < outputNames.size(); i++) {     remove(outputNames[i].c_str()); } m->Groups.clear(); return 0;  }
+               if (m->control_pressed) { outputTypes.clear(); delete input; for (int i = 0; i < outputNames.size(); i++) {     m->mothurRemove(outputNames[i]); } m->clearGroups(); return 0;  }
 
                //output error messages about any remaining user labels
                set<string>::iterator it;
@@ -367,7 +342,7 @@ int MatrixOutputCommand::execute(){
                        }
                }
                
-               if (m->control_pressed) { outputTypes.clear(); delete input;  for (int i = 0; i < outputNames.size(); i++) {    remove(outputNames[i].c_str()); } m->Groups.clear(); return 0;  }
+               if (m->control_pressed) { outputTypes.clear(); delete input;  for (int i = 0; i < outputNames.size(); i++) {    m->mothurRemove(outputNames[i]); } m->clearGroups(); return 0;  }
 
                //run last label if you need to
                if (needToRun == true)  {
@@ -379,10 +354,10 @@ int MatrixOutputCommand::execute(){
                        for (int i = 0; i < lookup.size(); i++) {  delete lookup[i];  } 
                }
                
-               if (m->control_pressed) { outputTypes.clear();  delete input;  for (int i = 0; i < outputNames.size(); i++) {   remove(outputNames[i].c_str()); } m->Groups.clear(); return 0;  }
+               if (m->control_pressed) { outputTypes.clear();  delete input;  for (int i = 0; i < outputNames.size(); i++) {   m->mothurRemove(outputNames[i]); } m->clearGroups(); return 0;  }
                
                //reset groups parameter
-               m->Groups.clear();  
+               m->clearGroups();  
                
                //set phylip file as new current phylipfile
                string current = "";
@@ -405,7 +380,7 @@ int MatrixOutputCommand::execute(){
        }
 }
 /***********************************************************/
-void MatrixOutputCommand::printSims(ostream& out) {
+void MatrixOutputCommand::printSims(ostream& out, vector< vector<float> >& simMatrix) {
        try {
                
                out.setf(ios::fixed, ios::floatfield); out.setf(ios::showpoint);
@@ -439,63 +414,216 @@ void MatrixOutputCommand::printSims(ostream& out) {
 /***********************************************************/
 int MatrixOutputCommand::process(vector<SharedRAbundVector*> thisLookup){
        try {
+               EstOutput data;
+               vector<SharedRAbundVector*> subset;
+               vector< vector<seqDist> > calcDists; calcDists.resize(matrixCalculators.size()); //one for each calc, this will be used to make .dist files
+               
        
-                               EstOutput data;
-                               vector<SharedRAbundVector*> subset;
-
-                               //for each calculator                                                                                           
-                               for(int i = 0 ; i < matrixCalculators.size(); i++) {
+               if(processors == 1){
+                       driver(thisLookup, 0, numGroups, calcDists);
+               }else{
+                       int process = 1;
+                       vector<int> processIDS;
+            
+                       #if defined (__APPLE__) || (__MACH__) || (linux) || (__linux)
+                       //loop through and create all the processes you want
+                       while (process != processors) {
+                               int pid = fork();
+                               
+                               if (pid > 0) {
+                                       processIDS.push_back(pid); 
+                                       process++;
+                               }else if (pid == 0){
+                                       driver(thisLookup, lines[process].start, lines[process].end, calcDists);   
                                        
-                                       //initialize simMatrix
-                                       simMatrix.clear();
-                                       simMatrix.resize(numGroups);
-                                       for (int p = 0; p < simMatrix.size(); p++)      {
-                                               for (int j = 0; j < simMatrix.size(); j++)      {
-                                                       simMatrix[p].push_back(0.0);
+                                       string tempdistFileName = m->getRootName(m->getSimpleName(sharedfile)) + toString(getpid()) + ".dist";
+                                       ofstream outtemp;
+                                       m->openOutputFile(tempdistFileName, outtemp);
+                                               
+                                       for (int i = 0; i < calcDists.size(); i++) {
+                                               outtemp << calcDists[i].size() << endl;
+                                                       
+                                               for (int j = 0; j < calcDists[i].size(); j++) {
+                                                       outtemp << calcDists[i][j].seq1 << '\t' << calcDists[i][j].seq2 << '\t' << calcDists[i][j].dist << 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);
+                               }
+                       }
+                       
+                       //parent do your part
+                       driver(thisLookup, lines[0].start, lines[0].end, calcDists);   
+                                               
+                       //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 tempdistFileName = m->getRootName(m->getSimpleName(sharedfile)) + toString(processIDS[i]) +  ".dist";
+                               ifstream intemp;
+                               m->openInputFile(tempdistFileName, intemp);
+                                       
+                               for (int k = 0; k < calcDists.size(); k++) {
+                                       int size = 0;
+                                       intemp >> size; m->gobble(intemp);
+                                               
+                                       for (int j = 0; j < size; j++) {
+                                               int seq1 = 0;
+                                               int seq2 = 0;
+                                               float dist = 1.0;
+                                                       
+                                               intemp >> seq1 >> seq2 >> dist;   m->gobble(intemp);
+                                                       
+                                               seqDist tempDist(seq1, seq2, dist);
+                                               calcDists[k].push_back(tempDist);
+                                       }
+                               }
+                               intemp.close();
+                               m->mothurRemove(tempdistFileName);
+                       }
+            #else
+            //////////////////////////////////////////////////////////////////////////////////////////////////////
+            //Windows version shared memory, so be careful when passing variables through the distSharedData 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<distSharedData*> 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;
+                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);
+                }
+                
+                //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.
+                distSharedData* tempSum = new distSharedData(m, lines[i].start, lines[i].end, Estimators, newLookup);
+                pDataArray.push_back(tempSum);
+                processIDS.push_back(i);
+                
+                hThreadArray[i-1] = CreateThread(NULL, 0, MyDistSharedThreadFunction, pDataArray[i-1], 0, &dwThreadIdArray[i-1]);   
+            }
+            
+            //parent do your part
+            driver(thisLookup, lines[0].start, lines[0].end, calcDists);   
+                       
+            //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 k = 0; k < calcDists.size(); k++) {
+                    int size = pDataArray[i]->calcDists[k].size();
+                    for (int j = 0; j < size; j++) {    calcDists[k].push_back(pDataArray[i]->calcDists[k][j]);    }
+                }
+                
+                CloseHandle(hThreadArray[i]);
+                delete pDataArray[i];
+            }
+
+            #endif
+               }
+
+               
+               
+               for (int i = 0; i < calcDists.size(); i++) {
+                       if (m->control_pressed) { break; }
+                               
+                       //initialize matrix
+                       vector< vector<float> > matrix; //square matrix to represent the distance
+                       matrix.resize(thisLookup.size());
+                       for (int k = 0; k < thisLookup.size(); k++) {  matrix[k].resize(thisLookup.size(), 0.0); }
+                               
+                       for (int j = 0; j < calcDists[i].size(); j++) {
+                               int row = calcDists[i][j].seq1;
+                               int column = calcDists[i][j].seq2;
+                               float dist = calcDists[i][j].dist;
+                                       
+                               matrix[row][column] = dist;
+                               matrix[column][row] = dist;
+                       }
+                       
+                       string distFileName = outputDir + m->getRootName(m->getSimpleName(sharedfile)) + matrixCalculators[i]->getName() + "." + thisLookup[0]->getLabel()  + "." + output + ".dist";
+                       outputNames.push_back(distFileName); outputTypes["phylip"].push_back(distFileName);
+                       ofstream outDist;
+                       m->openOutputFile(distFileName, outDist);
+                       outDist.setf(ios::fixed, ios::floatfield); outDist.setf(ios::showpoint);
+                       
+                       printSims(outDist, matrix);
+                       
+                       outDist.close();
+               }
+               
+               return 0;
+       }
+       catch(exception& e) {
+               m->errorOut(e, "MatrixOutputCommand", "process");
+               exit(1);
+       }
+}
+/**************************************************************************************************/
+int MatrixOutputCommand::driver(vector<SharedRAbundVector*> thisLookup, int start, int end, vector< vector<seqDist> >& calcDists) { 
+       try {
+               
+               vector<SharedRAbundVector*> subset;
+               for (int k = start; k < end; k++) { // pass cdd each set of groups to compare
+                       
+                       for (int l = 0; l < k; l++) {
                                
-                                       for (int k = 0; k < thisLookup.size(); k++) { 
-                                               for (int l = k; l < thisLookup.size(); l++) {
-                                                       if (k != l) { //we dont need to similiarity of a groups to itself
-                                                               //get estimated similarity between 2 groups
-                                                               
-                                                               if (m->control_pressed) { return 0; }
-                                                               
-                                                               subset.clear(); //clear out old pair of sharedrabunds
-                                                               //add new pair of sharedrabunds
-                                                               subset.push_back(thisLookup[k]); subset.push_back(thisLookup[l]); 
-                                                               
-                                                               //if this calc needs all groups to calculate the pair load all groups
-                                                               if (matrixCalculators[i]->getNeedsAll()) { 
-                                                                       //load subset with rest of lookup for those calcs that need everyone to calc for a pair
-                                                                       for (int w = 0; w < thisLookup.size(); w++) {
-                                                                               if ((w != k) && (w != l)) { subset.push_back(thisLookup[w]); }
-                                                                       }
-                                                               }
-                                                               
-                                                               data = matrixCalculators[i]->getValues(subset); //saves the calculator outputs
-                                                               //save values in similarity matrix
-                                                               simMatrix[k][l] = 1.0 - data[0];  //convert similiarity to distance
-                                                               simMatrix[l][k] = 1.0 - data[0];  //convert similiarity to distance
+                               if (k != l) { //we dont need to similiarity of a groups to itself
+                                       subset.clear(); //clear out old pair of sharedrabunds
+                                       //add new pair of sharedrabunds
+                                       subset.push_back(thisLookup[k]); subset.push_back(thisLookup[l]); 
+                                       
+                                       for(int i=0;i<matrixCalculators.size();i++) {
+                                               
+                                               //if this calc needs all groups to calculate the pair load all groups
+                                               if (matrixCalculators[i]->getNeedsAll()) { 
+                                                       //load subset with rest of lookup for those calcs that need everyone to calc for a pair
+                                                       for (int w = 0; w < thisLookup.size(); w++) {
+                                                               if ((w != k) && (w != l)) { subset.push_back(thisLookup[w]); }
                                                        }
                                                }
+                                               
+                                               vector<double> tempdata = matrixCalculators[i]->getValues(subset); //saves the calculator outputs
+                                               
+                                               if (m->control_pressed) { return 1; }
+                                               
+                                               seqDist temp(l, k, tempdata[0]);
+                                               calcDists[i].push_back(temp);
                                        }
-                                       
-                                       exportFileName = outputDir + m->getRootName(m->getSimpleName(sharedfile)) + matrixCalculators[i]->getName() + "." + thisLookup[0]->getLabel() + "." + output + ".dist";
-                                       m->openOutputFile(exportFileName, out);
-                                       outputNames.push_back(exportFileName); outputTypes["phylip"].push_back(exportFileName);
-                                       
-                                       printSims(out);
-                                       out.close();
-                                       
                                }
-
-                               return 0;
+                       }
+               }
                
+               return 0;
        }
        catch(exception& e) {
-               m->errorOut(e, "MatrixOutputCommand", "process");
+               m->errorOut(e, "MatrixOutputCommand", "driver");
                exit(1);
        }
 }