]> git.donarmstrong.com Git - mothur.git/blobdiff - parsimonycommand.cpp
precluster command finished
[mothur.git] / parsimonycommand.cpp
index 56305e2260d4af19050dfc0832add36c1208d13e..6ffb30b3caf18a695604c69c5556b2b7c5c2b50a 100644 (file)
 #include "parsimonycommand.h"
 
 /***********************************************************/
-ParsimonyCommand::ParsimonyCommand() {
+ParsimonyCommand::ParsimonyCommand(string option) {
        try {
                globaldata = GlobalData::getInstance();
+               abort = false;
+               Groups.clear();
                
-               //randomtree will tell us if user had their own treefile or if they just want the random distribution
-               randomtree = globaldata->getRandomTree();
+               //allow user to run help
+               if(option == "help") { help(); abort = true; }
                
-               //user has entered their own tree
-               if (randomtree == "") { 
-                       T = globaldata->gTree;
-                       tmap = globaldata->gTreemap;
-                       output = new ColumnFile(globaldata->getTreeFile()  +  ".parsimony");
-                       sumFile = globaldata->getTreeFile() + ".psummary";
-                       openOutputFile(sumFile, outSum);
-               }else { //user wants random distribution
-                       savetmap = globaldata->gTreemap;
-                       getUserInput();
-                       output = new ColumnFile(randomtree);
-               }
-               
-               //set users groups to analyze
-               util = new SharedUtil();
-               util->setGroups(globaldata->Groups, tmap->namesOfGroups, allGroups, numGroups, "unweighted");   //sets the groups the user wants to analyze
-               util->getCombos(groupComb, globaldata->Groups, numComp);
-               globaldata->setGroups("");
-               
-               if (numGroups == 1) { numComp++; groupComb.push_back(allGroups); }
+               else {
+                       //valid paramters for this command
+                       string Array[] =  {"random","groups","iters"};
+                       vector<string> myArray (Array, Array+(sizeof(Array)/sizeof(string)));
+                       
+                       OptionParser parser(option);
+                       map<string, string> parameters = parser.getParameters();
+                       
+                       ValidParameters validParameter;
                
-               convert(globaldata->getIters(), iters);  //how many random trees to generate
-               pars = new Parsimony(tmap);
-               counter = 0;
+                       //check to make sure all parameters are valid for command
+                       for (map<string,string>::iterator it = parameters.begin(); it != parameters.end(); it++) { 
+                               if (validParameter.isValidParameter(it->first, myArray, it->second) != true) {  abort = true;  }
+                       }
+                       
+                       randomtree = validParameter.validFile(parameters, "random", false);             if (randomtree == "not found") { randomtree = ""; }
+                       
+                       //are you trying to use parsimony without reading a tree or saying you want random distribution
+                       if (randomtree == "")  {
+                               if (globaldata->gTree.size() == 0) {
+                                       mothurOut("You must read a treefile and a groupfile or set the randomtree parameter to the output filename you wish, before you may execute the parsimony command."); mothurOutEndLine(); abort = true;  }
+                       }
+                                               
+                       //check for optional parameter and set defaults
+                       // ...at some point should added some additional type checking...
+                       groups = validParameter.validFile(parameters, "groups", false);                 
+                       if (groups == "not found") { groups = ""; globaldata->Groups.clear(); }
+                       else { 
+                               splitAtDash(groups, Groups);
+                               globaldata->Groups = Groups;
+                       }
+                               
+                       itersString = validParameter.validFile(parameters, "iters", false);                     if (itersString == "not found") { itersString = "1000"; }
+                       convert(itersString, iters); 
+                                               
+                       if (abort == false) {
+                               //randomtree will tell us if user had their own treefile or if they just want the random distribution
+                               //user has entered their own tree
+                               if (randomtree == "") { 
+                                       T = globaldata->gTree;
+                                       tmap = globaldata->gTreemap;
+                                       output = new ColumnFile(globaldata->getTreeFile()  +  ".parsimony", itersString);
+                                       sumFile = globaldata->getTreeFile() + ".psummary";
+                                       openOutputFile(sumFile, outSum);
+                               }else { //user wants random distribution
+                                       savetmap = globaldata->gTreemap;
+                                       getUserInput();
+                                       output = new ColumnFile(randomtree, itersString);
+                               }
+                               
+                               //set users groups to analyze
+                               util = new SharedUtil();
+                               util->setGroups(globaldata->Groups, tmap->namesOfGroups, allGroups, numGroups, "parsimony");    //sets the groups the user wants to analyze
+                               util->getCombos(groupComb, globaldata->Groups, numComp);
+                               
+                               if (numGroups == 1) { numComp++; groupComb.push_back(allGroups); }
+                               
+                               pars = new Parsimony(tmap);
+                               counter = 0;
+                               
+                       }
+                       
+               }
 
        }
        catch(exception& e) {
-               cout << "Standard Error: " << e.what() << " has occurred in the ParsimonyCommand class Function ParsimonyCommand. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
+               errorOut(e, "ParsimonyCommand", "ParsimonyCommand");
                exit(1);
        }
-       catch(...) {
-               cout << "An unknown error has occurred in the ParsimonyCommand class function ParsimonyCommand. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
+}
+
+//**********************************************************************************************************************
+
+void ParsimonyCommand::help(){
+       try {
+               mothurOut("The parsimony command can only be executed after a successful read.tree command, unless you use the random parameter.\n");
+               mothurOut("The parsimony command parameters are random, groups and iters.  No parameters are required.\n");
+               mothurOut("The groups parameter allows you to specify which of the groups in your groupfile you would like analyzed.  You must enter at least 1 valid group.\n");
+               mothurOut("The group names are separated by dashes.  The iters parameter allows you to specify how many random trees you would like compared to your tree.\n");
+               mothurOut("The parsimony command should be in the following format: parsimony(random=yourOutputFilename, groups=yourGroups, iters=yourIters).\n");
+               mothurOut("Example parsimony(random=out, iters=500).\n");
+               mothurOut("The default value for random is "" (meaning you want to use the trees in your inputfile, randomtree=out means you just want the random distribution of trees outputted to out.rd_parsimony),\n");
+               mothurOut("and iters is 1000.  The parsimony command output two files: .parsimony and .psummary their descriptions are in the manual.\n");
+               mothurOut("Note: No spaces between parameter labels (i.e. random), '=' and parameters (i.e.yourOutputFilename).\n\n");
+       }
+       catch(exception& e) {
+               errorOut(e, "ParsimonyCommand", "help");
                exit(1);
        }
 }
+
+
 /***********************************************************/
 int ParsimonyCommand::execute() {
        try {
+       
+               if (abort == true) { return 0; }
+       
                Progress* reading;
                reading = new Progress("Comparing to random:", iters);
                
@@ -77,7 +139,7 @@ int ParsimonyCommand::execute() {
                                for(int k = 0; k < numComp; k++) {
 
                                        //update uscoreFreq
-                                       it = uscoreFreq[k].find(userData[k]);
+                                       map<int,double>::iterator it = uscoreFreq[k].find(userData[k]);
                                        if (it == uscoreFreq[k].end()) {//new score
                                                uscoreFreq[k][userData[k]] = 1;
                                        }else{ uscoreFreq[k][userData[k]]++; }
@@ -103,8 +165,8 @@ int ParsimonyCommand::execute() {
                                        
                                for(int r = 0; r < numComp; r++) {
                                        //add trees pscore to map of scores
-                                       it2 = rscoreFreq[r].find(randomData[r]);
-                                       if (it2 != rscoreFreq[r].end()) {//already have that score
+                                       map<int,double>::iterator it = rscoreFreq[r].find(randomData[r]);
+                                       if (it != rscoreFreq[r].end()) {//already have that score
                                                rscoreFreq[r][randomData[r]]++;
                                        }else{//first time we have seen this score
                                                rscoreFreq[r][randomData[r]] = 1;
@@ -134,8 +196,8 @@ int ParsimonyCommand::execute() {
                        
                                for(int r = 0; r < numComp; r++) {
                                        //add trees pscore to map of scores
-                                       it2 = rscoreFreq[r].find(randomData[r]);
-                                       if (it2 != rscoreFreq[r].end()) {//already have that score
+                                       map<int,double>::iterator it = rscoreFreq[r].find(randomData[r]);
+                                       if (it != rscoreFreq[r].end()) {//already have that score
                                                rscoreFreq[r][randomData[r]]++;
                                        }else{//first time we have seen this score
                                                rscoreFreq[r][randomData[r]] = 1;
@@ -156,9 +218,9 @@ int ParsimonyCommand::execute() {
                        float rcumul = 0.0000;
                        float ucumul = 0.0000;
                        //this loop fills the cumulative maps and put 0.0000 in the score freq map to make it easier to print.
-                       for (it = validScores.begin(); it != validScores.end(); it++) { 
+                       for (map<int,double>::iterator it = validScores.begin(); it != validScores.end(); it++) { 
                                if (randomtree == "") {
-                                       it2 = uscoreFreq[a].find(it->first);
+                                       map<int,double>::iterator it2 = uscoreFreq[a].find(it->first);
                                        //user data has that score 
                                        if (it2 != uscoreFreq[a].end()) { uscoreFreq[a][it->first] /= T.size(); ucumul+= it2->second;  }
                                        else { uscoreFreq[a][it->first] = 0.0000; } //no user trees with that score
@@ -167,7 +229,7 @@ int ParsimonyCommand::execute() {
                                }
                        
                                //make rscoreFreq map and rCumul
-                               it2 = rscoreFreq[a].find(it->first);
+                               map<int,double>::iterator it2 = rscoreFreq[a].find(it->first);
                                //get percentage of random trees with that info
                                if (it2 != rscoreFreq[a].end()) {  rscoreFreq[a][it->first] /= iters; rcumul+= it2->second;  }
                                else { rscoreFreq[a][it->first] = 0.0000; } //no random trees with that score
@@ -195,8 +257,6 @@ int ParsimonyCommand::execute() {
                        globaldata->gTreemap = savetmap;
                }
                
-               //reset randomTree parameter to ""
-               globaldata->setRandomTree("");
                //reset groups parameter
                globaldata->Groups.clear(); 
                
@@ -204,11 +264,7 @@ int ParsimonyCommand::execute() {
                
        }
        catch(exception& e) {
-               cout << "Standard Error: " << e.what() << " has occurred in the ParsimonyCommand class Function execute. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
-               exit(1);
-       }
-       catch(...) {
-               cout << "An unknown error has occurred in the ParsimonyCommand class function execute. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
+               errorOut(e, "ParsimonyCommand", "execute");
                exit(1);
        }
 }
@@ -228,7 +284,7 @@ void ParsimonyCommand::printParsimonyFile() {
                for(int a = 0; a < numComp; a++) {
                        output->initFile(groupComb[a], tags);
                        //print each line
-                       for (it = validScores.begin(); it != validScores.end(); it++) { 
+                       for (map<int,double>::iterator it = validScores.begin(); it != validScores.end(); it++) { 
                                if (randomtree == "") {
                                        data.push_back(it->first);  data.push_back(uscoreFreq[a][it->first]); data.push_back(uCumul[a][it->first]); data.push_back(rscoreFreq[a][it->first]); data.push_back(rCumul[a][it->first]); 
                                }else{
@@ -241,11 +297,7 @@ void ParsimonyCommand::printParsimonyFile() {
                }
        }
        catch(exception& e) {
-               cout << "Standard Error: " << e.what() << " has occurred in the ParsimonyCommand class Function printParsimonyFile. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
-               exit(1);
-       }
-       catch(...) {
-               cout << "An unknown error has occurred in the ParsimonyCommand class function printParsimonyFile. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
+               errorOut(e, "ParsimonyCommand", "printParsimonyFile");
                exit(1);
        }
 }
@@ -254,7 +306,7 @@ void ParsimonyCommand::printUSummaryFile() {
        try {
                //column headers
                outSum << "Tree#" << '\t' << "Groups" << '\t'  <<  "ParsScore" << '\t' << "ParsSig" <<  endl;
-               cout << "Tree#" << '\t' << "Groups" << '\t'  <<  "ParsScore" << '\t' << "ParsSig" <<  endl;
+               mothurOut("Tree#\tGroups\tParsScore\tParsSig"); mothurOutEndLine();
                
                //format output
                outSum.setf(ios::fixed, ios::floatfield); outSum.setf(ios::showpoint);
@@ -264,11 +316,13 @@ void ParsimonyCommand::printUSummaryFile() {
                for (int i = 0; i< T.size(); i++) {
                        for(int a = 0; a < numComp; a++) {
                                if (UScoreSig[a][i] > (1/(float)iters)) {
-                                       outSum << setprecision(6) << i+1 << '\t' << groupComb[a]  << '\t' << userTreeScores[a][i] << setprecision(globaldata->getIters().length()) << '\t' << UScoreSig[a][i] << endl;
-                                       cout << setprecision(6) << i+1 << '\t' << groupComb[a]  << '\t' << userTreeScores[a][i] << setprecision(globaldata->getIters().length()) << '\t' << UScoreSig[a][i] << endl;
+                                       outSum << setprecision(6) << i+1 << '\t' << groupComb[a]  << '\t' << userTreeScores[a][i] << setprecision(itersString.length()) << '\t' << UScoreSig[a][i] << endl;
+                                       cout << setprecision(6) << i+1 << '\t' << groupComb[a]  << '\t' << userTreeScores[a][i] << setprecision(itersString.length()) << '\t' << UScoreSig[a][i] << endl;
+                                       mothurOutJustToLog(toString(i+1) + "\t" + groupComb[a] + "\t" + toString(userTreeScores[a][i]) + "\t" + toString(UScoreSig[a][i])); mothurOutEndLine();
                                }else {
-                                       outSum << setprecision(6) << i+1 << '\t' << groupComb[a] << '\t' << userTreeScores[a][i] << setprecision(globaldata->getIters().length())  << '\t' << "<" << (1/float(iters)) << endl;
-                                       cout << setprecision(6) << i+1 << '\t' << groupComb[a] << '\t' << userTreeScores[a][i] << setprecision(globaldata->getIters().length()) << '\t' << "<" << (1/float(iters)) << endl;
+                                       outSum << setprecision(6) << i+1 << '\t' << groupComb[a] << '\t' << userTreeScores[a][i] << setprecision(itersString.length())  << '\t' << "<" << (1/float(iters)) << endl;
+                                       cout << setprecision(6) << i+1 << '\t' << groupComb[a] << '\t' << userTreeScores[a][i] << setprecision(itersString.length()) << '\t' << "<" << (1/float(iters)) << endl;
+                                       mothurOutJustToLog(toString(i+1) + "\t" + groupComb[a] + "\t" + toString(userTreeScores[a][i]) + "\t" + toString((1/float(iters)))); mothurOutEndLine();
                                }
                        }
                }
@@ -276,11 +330,7 @@ void ParsimonyCommand::printUSummaryFile() {
                outSum.close();
        }
        catch(exception& e) {
-               cout << "Standard Error: " << e.what() << " has occurred in the ParsimonyCommand class Function printUSummaryFile. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
-               exit(1);
-       }
-       catch(...) {
-               cout << "An unknown error has occurred in the ParsimonyCommand class function printUSummaryFile. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
+               errorOut(e, "ParsimonyCommand", "printUSummaryFile");
                exit(1);
        }
 }
@@ -292,16 +342,18 @@ void ParsimonyCommand::getUserInput() {
                //create treemap
                tmap = new TreeMap();
 
-               cout << "Please enter the number of groups you would like to analyze: ";
+               mothurOut("Please enter the number of groups you would like to analyze: ");
                cin >> numGroups;
-                       
+               mothurOutJustToLog(toString(numGroups)); mothurOutEndLine();
+                               
                int num, count;
                count = 1;
                numEachGroup.resize(numGroups, 0);  
                
                for (int i = 1; i <= numGroups; i++) {
-                       cout << "Please enter the number of sequences in group " << i <<  ": ";
+                       mothurOut("Please enter the number of sequences in group " + toString(i) +  ": ");
                        cin >> num;
+                       mothurOutJustToLog(toString(num)); mothurOutEndLine();
                                
                        //set tmaps seqsPerGroup
                        tmap->seqsPerGroup[toString(i)] = num;
@@ -323,14 +375,11 @@ void ParsimonyCommand::getUserInput() {
                //memory leak prevention
                //if (globaldata->gTreemap != NULL) { delete globaldata->gTreemap;  }
                globaldata->gTreemap = tmap;
+               globaldata->Treenames = tmap->namesOfSeqs; 
                
        }
        catch(exception& e) {
-               cout << "Standard Error: " << e.what() << " has occurred in the ParsimonyCommand class Function getUserInput. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
-               exit(1);
-       }
-       catch(...) {
-               cout << "An unknown error has occurred in the ParsimonyCommand class function getUserInput. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
+               errorOut(e, "ParsimonyCommand", "getUserInput");
                exit(1);
        }
 }