]> git.donarmstrong.com Git - mothur.git/blobdiff - parsimonycommand.cpp
put back no command and worked on libshuff
[mothur.git] / parsimonycommand.cpp
index e741b5e914770e34b49ac5c244e31b041f8520d2..493af9b867e5542ee4c4c93acfa6c40a9c95fc5a 100644 (file)
@@ -22,20 +22,21 @@ ParsimonyCommand::ParsimonyCommand() {
                        T = globaldata->gTree;
                        tmap = globaldata->gTreemap;
                        parsFile = globaldata->getTreeFile() + ".parsimony";
-                       openOutputFile(parsFile, out);
+                       parsFileout = globaldata->getTreeFile() + "temp" + ".parsimony";
                        sumFile = globaldata->getTreeFile() + ".psummary";
                        openOutputFile(sumFile, outSum);
-                       //set users groups to analyze
-                       setGroups();
                }else { //user wants random distribution
                        savetmap = globaldata->gTreemap;
                        getUserInput();
-                       parsFile = randomtree + ".rd_parsimony";
-                       openOutputFile(parsFile, out);
+                       parsFile = randomtree;
+                       parsFileout = globaldata->getTreeFile() + "temp";
                }
                
+               //set users groups to analyze
+               setGroups();
                convert(globaldata->getIters(), iters);  //how many random trees to generate
                pars = new Parsimony(tmap);
+               counter = 0;
 
        }
        catch(exception& e) {
@@ -50,28 +51,24 @@ ParsimonyCommand::ParsimonyCommand() {
 /***********************************************************/
 int ParsimonyCommand::execute() {
        try {
-               
+       
+               //get pscore for users tree
+               userData.resize(numComp,0);  //data = AB, AC, BC, ABC.
+               randomData.resize(numComp,0);  //data = AB, AC, BC, ABC.
+               rscoreFreq.resize(numComp);  
+               uscoreFreq.resize(numComp);  
+               rCumul.resize(numComp);  
+               uCumul.resize(numComp);  
+               userTreeScores.resize(numComp);  
+               UScoreSig.resize(numComp); 
                                
                if (randomtree == "") {
-                       //get pscore for users tree
-                       userData.resize(numComp,0);  //data = AB, AC, BC, ABC.
-                       randomData.resize(numComp,0);  //data = AB, AC, BC, ABC.
-                       rscoreFreq.resize(numComp);  
-                       uscoreFreq.resize(numComp);  
-                       rCumul.resize(numComp);  
-                       uCumul.resize(numComp);  
-                       validScores.resize(numComp); 
-                       userTreeScores.resize(numComp);  
-                       UScoreSig.resize(numComp); 
-
                        //get pscores for users trees
                        for (int i = 0; i < T.size(); i++) {
-                               cout << "Processing tree " << i+1 << endl;
                                userData = pars->getValues(T[i]);  //data = AB, AC, BC, ABC.
                                
                                //output scores for each combination
                                for(int k = 0; k < numComp; k++) {
-                                       cout << "Tree " << i+1 << " Combination " << groupComb[k] << " parsimony score = " << userData[k] << endl;
                                        //update uscoreFreq
                                        it = uscoreFreq[k].find(userData[k]);
                                        if (it == uscoreFreq[k].end()) {//new score
@@ -79,7 +76,7 @@ int ParsimonyCommand::execute() {
                                        }else{ uscoreFreq[k][userData[k]]++; }
                                        
                                        //add users score to valid scores
-                                       validScores[k][userData[k]] = userData[k];
+                                       validScores[userData[k]] = userData[k];
                                        
                                        //save score for summary file
                                        userTreeScores[k].push_back(userData[k]);
@@ -105,7 +102,7 @@ int ParsimonyCommand::execute() {
                                        }
                        
                                        //add randoms score to validscores
-                                       validScores[r][randomData[r]] = randomData[r];
+                                       validScores[randomData[r]] = randomData[r];
                                }
                                
                                delete randT;
@@ -130,26 +127,25 @@ int ParsimonyCommand::execute() {
                                        }
                        
                                        //add randoms score to validscores
-                                       validScores[r][randomData[r]] = randomData[r];
+                                       validScores[randomData[r]] = randomData[r];
                                }
                                
                                delete randT;
                        }
                }
                
-               float rcumul = 0.0000;
-               float ucumul = 0.0000;
-               
                for(int a = 0; a < numComp; a++) {
-               //this loop fills the cumulative maps and put 0.0000 in the score freq map to make it easier to print.
-                       for (it = validScores[a].begin(); it != validScores[a].end(); it++) { 
+                       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++) { 
                                if (randomtree == "") {
                                        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
                                        //make uCumul map
-                                       uCumul[a][it->first] = ucumul-a;
+                                       uCumul[a][it->first] = ucumul;
                                }
                        
                                //make rscoreFreq map and rCumul
@@ -157,7 +153,7 @@ int ParsimonyCommand::execute() {
                                //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
-                               rCumul[a][it->first] = rcumul-a;
+                               rCumul[a][it->first] = rcumul;
                        }
                        
                        //find the signifigance of each user trees score when compared to the random trees and save for printing the summary file
@@ -167,10 +163,14 @@ int ParsimonyCommand::execute() {
                }
                
                printParsimonyFile();
-               printUSummaryFile();
+               if (randomtree == "") { printUSummaryFile(); }
                
                //reset globaldata's treemap if you just did random distrib
-               if (randomtree != "") { globaldata->gTreemap = savetmap; }
+               if (randomtree != "") {
+                       //memory leak prevention
+                       //if (globaldata->gTreemap != NULL) { delete globaldata->gTreemap;  }
+                       globaldata->gTreemap = savetmap;
+               }
                
                //reset randomTree parameter to ""
                globaldata->setRandomTree("");
@@ -193,28 +193,29 @@ int ParsimonyCommand::execute() {
 /***********************************************************/
 void ParsimonyCommand::printParsimonyFile() {
        try {
-               //column headers
-               if (randomtree == "") {
-                       out << "Comb" << '\t' << "Score" << '\t' << "UserFreq" << '\t' << "UserCumul" << '\t' << "RandFreq" << '\t' << "RandCumul" << endl;
-               }else {
-                       out << "Comb" << '\t' << "Score" << '\t' << "RandFreq" << '\t' << "RandCumul" << endl;
-               }
+               vector<double> data;
                
                //format output
                out.setf(ios::fixed, ios::floatfield); out.setf(ios::showpoint);
-               
+
                for(int a = 0; a < numComp; a++) {
+                       initFile(groupComb[a]);
                        //print each line
-                       for (it = validScores[a].begin(); it != validScores[a].end(); it++) { 
+                       for (it = validScores.begin(); it != validScores.end(); it++) { 
                                if (randomtree == "") {
-                                       out << setprecision(6)  << groupComb[a] << '\t' << it->first << '\t' << '\t'<< uscoreFreq[a][it->first] << '\t' << uCumul[a][it->first] << '\t' << rscoreFreq[a][it->first] << '\t' << rCumul[a][it->first] << endl
+                                       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{
-                                       out << setprecision(6) << groupComb[a] << '\t' << it->first << '\t' << '\t' << rscoreFreq[a][it->first] << '\t' << rCumul[a][it->first] << endl;        
+                                       data.push_back(it->first);  data.push_back(rscoreFreq[a][it->first]); data.push_back(rCumul[a][it->first]); 
                                }
+                               output(data);
+                               data.clear();
                        } 
+                       resetFile();
                }
-               out.close();
                
+               out.close();
+               inFile.close();
+               remove(parsFileout.c_str());
        }
        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";
@@ -229,7 +230,8 @@ void ParsimonyCommand::printParsimonyFile() {
 void ParsimonyCommand::printUSummaryFile() {
        try {
                //column headers
-               outSum << "Tree#" << '\t' << "Comb" << '\t'  <<  "ParsScore" << '\t' << '\t' << "ParsSig" <<  endl;
+               outSum << "Tree#" << '\t' << "Groups" << '\t'  <<  "ParsScore" << '\t' << "ParsSig" <<  endl;
+               cout << "Tree#" << '\t' << "Groups" << '\t'  <<  "ParsScore" << '\t' << "ParsSig" <<  endl;
                
                //format output
                outSum.setf(ios::fixed, ios::floatfield); outSum.setf(ios::showpoint);
@@ -238,7 +240,13 @@ void ParsimonyCommand::printUSummaryFile() {
                //print each line
                for (int i = 0; i< T.size(); i++) {
                        for(int a = 0; a < numComp; a++) {
-                               outSum << setprecision(6) << i+1 << '\t' << groupComb[a] << '\t' << '\t' << userTreeScores[a][i] << '\t' << UScoreSig[a][i] << endl;
+                               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;
+                               }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;
+                               }
                        }
                }
                
@@ -274,6 +282,7 @@ void ParsimonyCommand::getUserInput() {
                                
                        //set tmaps seqsPerGroup
                        tmap->seqsPerGroup[toString(i)] = num;
+                       tmap->namesOfGroups.push_back(toString(i));
                        
                        //set tmaps namesOfSeqs
                        for (int j = 0; j < num; j++) {
@@ -288,6 +297,8 @@ void ParsimonyCommand::getUserInput() {
                getline(cin, s);
                
                //save tmap for later
+               //memory leak prevention
+               //if (globaldata->gTreemap != NULL) { delete globaldata->gTreemap;  }
                globaldata->gTreemap = tmap;
                
        }
@@ -324,26 +335,31 @@ void ParsimonyCommand::setGroups() {
                                        for (int i = 0; i < tmap->namesOfGroups.size(); i++) {
                                                globaldata->Groups.push_back(tmap->namesOfGroups[i]);
                                                numGroups++;
-                                               allGroups += tmap->namesOfGroups[i];
+                                               allGroups += tmap->namesOfGroups[i] + "-";
                                        }
+                                       allGroups = allGroups.substr(0, allGroups.length()-1);
                                }else {
                                        for (int i = 0; i < globaldata->Groups.size(); i++) {
-                                               allGroups += tmap->namesOfGroups[i];
+                                               allGroups += globaldata->Groups[i] + "-";
                                                numGroups++;
                                        }
+                                       allGroups = allGroups.substr(0, allGroups.length()-1);
                                }
                        }else{//user has enter "all" and wants the default groups
+                               globaldata->Groups.clear();
                                for (int i = 0; i < tmap->namesOfGroups.size(); i++) {
                                        globaldata->Groups.push_back(tmap->namesOfGroups[i]);
                                        numGroups++;
-                                       allGroups += tmap->namesOfGroups[i];
+                                       allGroups += tmap->namesOfGroups[i] + "-";
                                }
+                               allGroups = allGroups.substr(0, allGroups.length()-1);
                                globaldata->setGroups("");
                        }
                }else {
                        for (int i = 0; i < tmap->namesOfGroups.size(); i++) {
-                               allGroups += tmap->namesOfGroups[i];
+                               allGroups += tmap->namesOfGroups[i] + "-";
                        }
+                       allGroups = allGroups.substr(0, allGroups.length()-1);
                        numGroups = 1;
                }
                
@@ -351,14 +367,16 @@ void ParsimonyCommand::setGroups() {
                numComp = 0;
                for (int r=0; r<numGroups; r++) { 
                        for (int l = r+1; l < numGroups; l++) {
-                               groupComb.push_back(globaldata->Groups[r]+globaldata->Groups[l]);
+                               groupComb.push_back(globaldata->Groups[r]+ "-" +globaldata->Groups[l]);
                                numComp++;
                        }
                }
                
                //ABC
-               groupComb.push_back(allGroups);
-               numComp++;
+               if (numComp != 1) {
+                       groupComb.push_back(allGroups);
+                       numComp++;
+               }
                
        }
        catch(exception& e) {
@@ -373,4 +391,100 @@ void ParsimonyCommand::setGroups() {
 }
 /*****************************************************************/
 
+void ParsimonyCommand::initFile(string label){
+       try {
+               if(counter != 0){
+                       openOutputFile(parsFileout, out);
+                       openInputFile(parsFile, inFile);
+
+                       string inputBuffer;
+                       getline(inFile, inputBuffer);
+                       
+                       if (randomtree == "") {
+                               out <<  inputBuffer << '\t' << label + "Score" << '\t' << label + "UserFreq" << '\t' << label + "UserCumul" << '\t' << label + "RandFreq" << '\t' << label + "RandCumul" << endl;
+                       }else {
+                               out <<  inputBuffer << '\t' << "Score" << '\t' << "RandFreq" << '\t' << "RandCumul" << endl;
+                       }
+               }else{
+                       openOutputFile(parsFileout, out);
+                       //column headers
+                       if (randomtree == "") {
+                               out << label + "Score" << '\t' << label + "UserFreq" << '\t' << label + "UserCumul" << '\t' << label + "RandFreq" << '\t' << label + "RandCumul" << endl;
+                       }else {
+                               out << "Score" << '\t' << "RandFreq" << '\t' << "RandCumul" << endl;
+                       }
+               }
+
+               out.setf(ios::fixed, ios::floatfield);
+               out.setf(ios::showpoint);
+       }
+       catch(exception& e) {
+               cout << "Standard Error: " << e.what() << " has occurred in the ParsimonyCommand class Function initFile. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
+               exit(1);
+       }
+       catch(...) {
+               cout << "An unknown error has occurred in the ParsimonyCommand class function initFile. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
+               exit(1);
+       }
+}
+
+/***********************************************************************/
+
+void ParsimonyCommand::output(vector<double> data){
+       try {
+               if(counter != 0){               
+                       string inputBuffer;
+                       getline(inFile, inputBuffer);
+               
+                       if (randomtree == "") {
+                               out << inputBuffer << '\t' << setprecision(6) << data[0] << setprecision(globaldata->getIters().length())  << '\t' << data[1] << '\t' << data[2] << '\t' << data[3] << '\t' << data[4] << endl;
+                       }else{
+                               out << inputBuffer << '\t' << setprecision(6) << data[0] << setprecision(globaldata->getIters().length())  << '\t' << data[1] << '\t' << data[2] << endl;
+                       }
+               }
+               else{
+                       if (randomtree == "") {
+                               out << setprecision(6) << data[0] << setprecision(globaldata->getIters().length())  << '\t' << data[1] << '\t' << data[2] << '\t' << data[3] << '\t' << data[4] << endl;
+                       }else{
+                               out << setprecision(6) << data[0] << setprecision(globaldata->getIters().length())  << '\t' << data[1] << '\t' << data[2] << endl;
+                       }
+               }
+
+       }
+       catch(exception& e) {
+               cout << "Standard Error: " << e.what() << " has occurred in the ParsimonyCommand class Function output. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
+               exit(1);
+       }
+       catch(...) {
+               cout << "An unknown error has occurred in the ParsimonyCommand class function output. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
+               exit(1);
+       }
+}
+
+/***********************************************************************/
+
+void ParsimonyCommand::resetFile(){
+       try {
+               if(counter != 0){
+                       out.close();
+                       inFile.close();
+               }
+               else{
+                       out.close();
+               }
+               counter = 1;
+               
+               remove(parsFile.c_str());
+               rename(parsFileout.c_str(), parsFile.c_str());
+       }
+       catch(exception& e) {
+               cout << "Standard Error: " << e.what() << " has occurred in the ParsimonyCommand class Function resetFile. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
+               exit(1);
+       }
+       catch(...) {
+               cout << "An unknown error has occurred in the ParsimonyCommand class function resetFile. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
+               exit(1);
+       }       
+}
+