]> git.donarmstrong.com Git - mothur.git/commitdiff
fixed unweighted calculator
authorwestcott <westcott>
Fri, 20 Feb 2009 17:34:55 +0000 (17:34 +0000)
committerwestcott <westcott>
Fri, 20 Feb 2009 17:34:55 +0000 (17:34 +0000)
14 files changed:
errorchecking.cpp
errorchecking.h
parsimony.cpp
parsimony.h
parsimonycommand.cpp
parsimonycommand.h
tree.cpp
tree.h
unifracunweightedcommand.cpp
unifracunweightedcommand.h
unweighted.cpp
unweighted.h
utilities.hpp
validparameter.cpp

index dccb679cd6638d6de391154b328301e4f7c17df9..ff94c8f4f765a8e0041604416e3b6807b423c05a 100644 (file)
@@ -36,9 +36,12 @@ void ErrorCheck::refresh() {
        cutoff = globaldata->getCutOff();
        format = globaldata->getFormat();
        method = globaldata->getMethod();
+       randomtree = globaldata->getRandomTree();
+       sharedfile = globaldata->getSharedFile();
+
 
        
-       string p[] = {
+/*     string p[] = {
                "phylip",              //0
                "column",             //1
                "list",               //2
@@ -155,10 +158,8 @@ void ErrorCheck::refresh() {
        intParams[p[13]] = ipv2;
        intParams[p[14]] = ipv3;
        intParams[p[17]] = ipv4;
-       intParams[p[26]] = ipv5;
+       intParams[p[26]] = ipv5;   */
        
-       randomtree = globaldata->getRandomTree();
-       sharedfile = globaldata->getSharedFile();
 }
 
 /*******************************************************/
@@ -203,15 +204,15 @@ bool ErrorCheck::checkInput(string input) {
                                
                                //is it a valid parameter
                                if (validParameter->isValidParameter(parameter) != true) { return false; }
-                               if(!validCommandParameter(parameter,commandName)) { 
-                                       cout << "'" << parameter << "' is not a valid parameter for the " << commandName << " command.\n";
-                                       return false; 
-                               }
-                               if(!validParameterValue(value, parameter)) {
-                                       if(parameter.compare("precision") == 0)
-                                               cout << "The precision parameter can only take powers of 10 as a value (e.g. 10,1000,1000, etc.)\n";
-                                       else {
-                                       vector<double> bounds = intParams[parameter];
+                               //if(!validCommandParameter(parameter,commandName)) { 
+                               //      cout << "'" << parameter << "' is not a valid parameter for the " << commandName << " command.\n";
+                               //      return false; 
+                               //}
+                               //if(!validParameterValue(value, parameter)) {
+                               //      if(parameter.compare("precision") == 0)
+                               //              cout << "The precision parameter can only take powers of 10 as a value (e.g. 10,1000,1000, etc.)\n";
+                               //      else {
+                               /*      vector<double> bounds = intParams[parameter];
                                        double a = bounds.at(0);
                                        double b = bounds.at(1);
                                        double c = bounds.at(2);
@@ -243,7 +244,7 @@ bool ErrorCheck::checkInput(string input) {
                                        }
                                        }
                                        return false;
-                               }
+                               } */
 
                                if (parameter == "phylip" )             { phylipfile = value; }
                                if (parameter == "column" )             { columnfile = value; }
@@ -275,11 +276,11 @@ bool ErrorCheck::checkInput(string input) {
                                splitAtEquals(parameter, value);
                                //is it a valid parameter
                                if (validParameter->isValidParameter(parameter) != true) { return false; }
-                               if(!validCommandParameter(parameter,commandName)) { 
-                                       cout << "'" << parameter << "' is not a valid parameter for the " << commandName << " command.\n";
-                                       return false; 
-                               }
-                               if(!validParameterValue(value, parameter)) {
+                       //      if(!validCommandParameter(parameter,commandName)) { 
+                       //              cout << "'" << parameter << "' is not a valid parameter for the " << commandName << " command.\n";
+                       //              return false; 
+                       //      }
+                       /*      if(!validParameterValue(value, parameter)) {
                                        if(parameter.compare("precision") == 0)
                                                cout << "The precision parameter can only take powers of 10 as a value (e.g. 10,1000,1000, etc.)\n";
                                        else {
@@ -315,7 +316,7 @@ bool ErrorCheck::checkInput(string input) {
                                        }
                                        }
                                        return false;
-                               }
+                               }*/
                                if (parameter == "phylip" )             { phylipfile = value; }
                                if (parameter == "column" )             { columnfile = value; }                         
                                if (parameter == "list" )               { listfile = value; }
@@ -485,7 +486,7 @@ void ErrorCheck::validateReadFiles() {
 }
 /*******************************************************/
 
-/******************************************************/
+/******************************************************
 //This function checks to see if the given paramter
 //is a valid paramter for the given command.
 bool ErrorCheck::validCommandParameter(string parameter, string commandName) {
@@ -507,7 +508,7 @@ bool ErrorCheck::validCommandParameter(string parameter, string commandName) {
 }
 /*******************************************************/
 
-/******************************************************/
+/******************************************************
 //This function checks to see if the given paramter value
 //is convertable into an int if that parameter requires it.
 bool ErrorCheck::validParameterValue(string value, string parameter) {
index 8bf68578a44330ce9c3c0c595c0200486b5b41c4..50cc00511a08fa71754d2dbf51fd76406d2925ea 100644 (file)
@@ -26,8 +26,8 @@ class ErrorCheck {
                ValidCommands* validCommand;
                ValidParameters* validParameter;
                void validateReadFiles();
-               bool validCommandParameter(string, string);
-               bool validParameterValue(string, string);
+       //      bool validCommandParameter(string, string);
+       //      bool validParameterValue(string, string);
                void validateReadDist();
                void validateReadPhil();
                void validateParseFiles();
@@ -40,9 +40,9 @@ class ErrorCheck {
                bool errorFree;
 
                vector<string> sharedGroups;
-               map <string, vector<string> > commandParameters;
-               map <string, vector<double> > intParams;
-               double piSent;
+       //      map <string, vector<string> > commandParameters;
+       ///     map <string, vector<double> > intParams;
+       //      double piSent;
 
 };
 #endif
index 2a673b0cb622bd97b8f91dc4f1c2e9c69de19e8d..6239f75f08a59c9a5f690401e15ef93c19fe8f50 100644 (file)
@@ -28,7 +28,7 @@ EstOutput Parsimony::getValues(Tree* t) {
 
                        //add in all the groups the users wanted
                        for (it = t->tree[i].pGroups.begin(); it != t->tree[i].pGroups.end(); it++) {
-                               if (inUsersGroups(it->first) == true) {  iSize++;  }
+                               if (inUsersGroups(it->first, globaldata->Groups) == true) {  iSize++;  }
                        }
 
                        //if that leaves no groups give it 1 so it will cause no change to parent
@@ -37,7 +37,7 @@ EstOutput Parsimony::getValues(Tree* t) {
                        //add in all the groups the users wanted
                        for (it = t->tree[rc].pGroups.begin(); it != t->tree[rc].pGroups.end(); it++) {
 
-                               if (inUsersGroups(it->first) == true) {  rcSize++;  }
+                               if (inUsersGroups(it->first, globaldata->Groups) == true) {  rcSize++;  }
                        }
                        
                        //if that leaves no groups give it 1 so it will cause no change to parent
@@ -47,7 +47,7 @@ EstOutput Parsimony::getValues(Tree* t) {
                        //add in all the groups the users wanted
                        for (it = t->tree[lc].pGroups.begin(); it != t->tree[lc].pGroups.end(); it++) {
 
-                               if (inUsersGroups(it->first) == true) {  lcSize++;  }
+                               if (inUsersGroups(it->first, globaldata->Groups) == true) {  lcSize++;  }
                        }
                        
                        //if that leaves no groups give it 1 so it will cause no change to parent
@@ -75,22 +75,4 @@ EstOutput Parsimony::getValues(Tree* t) {
        }
 
 }
-/**************************************************************************************************/
 
-bool Parsimony::inUsersGroups(string groupname) {
-       try {
-               for (int i = 0; i < globaldata->Groups.size(); i++) {
-                       if (groupname == globaldata->Groups[i]) { return true; }
-               }
-               return false;
-       }
-       catch(exception& e) {
-               cout << "Standard Error: " << e.what() << " has occurred in the Parsimony class Function inUsersGroups. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
-               exit(1);
-       }
-       catch(...) {
-               cout << "An unknown error has occurred in the Parsimony class function inUsersGroups. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
-               exit(1);
-       }
-}
-/**************************************************************************************************/
index 8955ef7a2253ec5a1f2985d09198431a66a8db2c..ebe3d8d89d4ecf2174dc84f8796245a07f4384c9 100644 (file)
@@ -28,7 +28,6 @@ class Parsimony : public TreeCalculator  {
                GlobalData* globaldata;
                EstOutput data;
                TreeMap* tmap;
-               bool inUsersGroups(string);
                map<string, int>::iterator it;
 };
 
index fe5a1dd251899818715cb52c14ac7c73a40630b2..554d24bd93f76c94a223d33ce24159babd0f3e78 100644 (file)
@@ -27,31 +27,8 @@ ParsimonyCommand::ParsimonyCommand() {
                        openOutputFile(sumFile, outSum);
                        distFile = globaldata->getTreeFile() + ".pdistrib";
                        openOutputFile(distFile, outDist);
-                       
-                       //if the user has not entered specific groups to analyze then do them all
-                       if (globaldata->Groups.size() != 0) {
-                               //check that groups are valid
-                               for (int i = 0; i < globaldata->Groups.size(); i++) {
-                                       if (tmap->isValidGroup(globaldata->Groups[i]) != true) {
-                                               cout << globaldata->Groups[i] << " is not a valid group, and will be disregarded." << endl;
-                                               // erase the invalid group from globaldata->Groups
-                                               globaldata->Groups.erase (globaldata->Groups.begin()+i);
-                                       }
-                               }
-                       
-                               //if the user only entered invalid groups
-                               if (globaldata->Groups.size() == 0) { 
-                                       cout << "When using the groups parameter you must have at least 1 valid group. I will run the command using all the groups in your groupfile." << endl; 
-                                       for (int i = 0; i < tmap->namesOfGroups.size(); i++) {
-                                               globaldata->Groups.push_back(tmap->namesOfGroups[i]);
-                                       }
-                               }               
-                       }else {
-                               for (int i = 0; i < tmap->namesOfGroups.size(); i++) {
-                                       globaldata->Groups.push_back(tmap->namesOfGroups[i]);
-                               }
-                       }
-
+                       //set users groups to analyze
+                       setGroups();
                }else { //user wants random distribution
                        savetmap = globaldata->gTreemap;
                        getUserInput();
@@ -313,3 +290,43 @@ void ParsimonyCommand::getUserInput() {
 }
 /***********************************************************/
 
+void ParsimonyCommand::setGroups() {
+       try {
+               //if the user has not entered specific groups to analyze then do them all
+               if (globaldata->Groups.size() != 0) {
+                       //check that groups are valid
+                       for (int i = 0; i < globaldata->Groups.size(); i++) {
+                               if (tmap->isValidGroup(globaldata->Groups[i]) != true) {
+                                       cout << globaldata->Groups[i] << " is not a valid group, and will be disregarded." << endl;
+                                       // erase the invalid group from globaldata->Groups
+                                       globaldata->Groups.erase (globaldata->Groups.begin()+i);
+                               }
+                       }
+                       
+                       //if the user only entered invalid groups
+                       if (globaldata->Groups.size() == 0) { 
+                               cout << "When using the groups parameter you must have at least 1 valid group. I will run the command using all the groups in your groupfile." << endl; 
+                               for (int i = 0; i < tmap->namesOfGroups.size(); i++) {
+                                       globaldata->Groups.push_back(tmap->namesOfGroups[i]);
+                               }
+                       }
+                                       
+               }else {
+                       for (int i = 0; i < tmap->namesOfGroups.size(); i++) {
+                               globaldata->Groups.push_back(tmap->namesOfGroups[i]);
+                       }
+               }
+       }
+       catch(exception& e) {
+               cout << "Standard Error: " << e.what() << " has occurred in the ParsimonyCommand class Function setGroups. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
+               exit(1);
+       }
+       catch(...) {
+               cout << "An unknown error has occurred in the ParsimonyCommand class function setGroups. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
+               exit(1);
+       }               
+
+}
+/*****************************************************************/
+
+
index bd52fa6a74b66a0700316ebb069519582d0b8733..326c0c9170ceae07e7061d31c46a25653f52cb7d 100644 (file)
@@ -51,6 +51,7 @@ class ParsimonyCommand : public Command {
                void printParsimonyFile();  
                void printUSummaryFile();
                void getUserInput();
+               void setGroups();
                
 };
 
index 94ab326d8f6ba32ededf4f5b2240dc5852e4816c..2ad46f044f71556a9726c93a96d4fcd5ae42a6b6 100644 (file)
--- a/tree.cpp
+++ b/tree.cpp
@@ -231,33 +231,51 @@ map<string,int> Tree::mergeGcounts(int position) {
 
 void Tree::randomLabels() {
        try {
+               
+               //set up the groups the user wants to include
+               setGroups();
+               
                for(int i=numLeaves-1;i>=0;i--){
                        if(tree[i].pGroups.size() == 0){
                                continue;
                        }
-               
+                       
                        int escape = 1;
                        int z;
                
                        while(escape == 1){
+                               //get random index to switch with
                                z = int((float)(i+1) * (float)(rand()) / ((float)RAND_MAX+1.0));        
                        
                                if(tree[z].pGroups.size() != 0){
                                        escape = 0;
                                }
                        }
+                       
+                       //you only want to randomize the nodes that are from a group the user wants analyzed, so
+                       //if either of the leaf nodes you are about to switch are not in the users groups then you don't want to switch them.
+                       bool treez, treei;
+                       
+                       //leaves have only one group so you can just set it to begin()
+                       it = tree[z].pGroups.begin();
+                       treez = inUsersGroups(it->first, globaldata->Groups);
+                       
+                       it = tree[i].pGroups.begin();
+                       treei = inUsersGroups(it->first, globaldata->Groups);
+                       
+                       if ((treez == true) && (treei == true)) {
+                               //switches node i and node z's info.
+                               map<string,int> lib_hold = tree[z].pGroups;
+                               tree[z].pGroups = (tree[i].pGroups);
+                               tree[i].pGroups = (lib_hold);
                
+                               tree[z].setGroup(tree[z].pGroups.begin()->first);
+                               tree[i].setGroup(tree[i].pGroups.begin()->first);
                
-                       map<string,int> lib_hold = tree[z].pGroups;
-                       tree[z].pGroups = (tree[i].pGroups);
-                       tree[i].pGroups = (lib_hold);
-               
-                       tree[z].setGroup(tree[z].pGroups.begin()->first);
-                       tree[i].setGroup(tree[i].pGroups.begin()->first);
-               
-                       map<string,int> gcount_hold = tree[z].pcount;
-                       tree[z].pcount = (tree[i].pcount);
-                       tree[i].pcount = (gcount_hold);
+                               map<string,int> gcount_hold = tree[z].pcount;
+                               tree[z].pcount = (tree[i].pcount);
+                               tree[i].pcount = (gcount_hold);
+                       }
                }
        }
        catch(exception& e) {
@@ -416,6 +434,41 @@ void Tree::printBranch(int node) {
 }
 
 /*****************************************************************/
+void Tree::setGroups() {
+       try {
+               //if the user has not entered specific groups to analyze then do them all
+               if (globaldata->Groups.size() != 0) {
+                       //check that groups are valid
+                       for (int i = 0; i < globaldata->Groups.size(); i++) {
+                               if (globaldata->gTreemap->isValidGroup(globaldata->Groups[i]) != true) {
+                                       cout << globaldata->Groups[i] << " is not a valid group, and will be disregarded." << endl;
+                                       // erase the invalid group from globaldata->Groups
+                                       globaldata->Groups.erase (globaldata->Groups.begin()+i);
+                               }
+                       }
+                       
+                       //if the user only entered invalid groups
+                       if (globaldata->Groups.size() == 0) { 
+                               cout << "When using the groups parameter you must have at least 1 valid group. I will run the command using all the groups in your groupfile." << endl; 
+                               for (int i = 0; i < globaldata->gTreemap->namesOfGroups.size(); i++) {
+                                       globaldata->Groups.push_back(globaldata->gTreemap->namesOfGroups[i]);
+                               }
+                       }
+                                       
+               }else {
+                       for (int i = 0; i < globaldata->gTreemap->namesOfGroups.size(); i++) {
+                               globaldata->Groups.push_back(globaldata->gTreemap->namesOfGroups[i]);
+                       }
+               }
+       }
+       catch(exception& e) {
+               cout << "Standard Error: " << e.what() << " has occurred in the Tree class Function setGroups. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
+               exit(1);
+       }
+       catch(...) {
+               cout << "An unknown error has occurred in the Tree class function setGroups. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
+               exit(1);
+       }               
 
-
+}
 
diff --git a/tree.h b/tree.h
index d7dcc3ddd790e0894347aa1bcb37d2da9608a567..7157a6f0e43993e8278591d3422caee1eeaa8728 100644 (file)
--- a/tree.h
+++ b/tree.h
@@ -51,6 +51,7 @@ class Tree {
                void randomLabels();
                int findRoot();  //return index of root node
                void printBranch(int);  //recursively print out tree
+               void setGroups();
 };
 
 #endif
index bf2129c8b2bbc2e8e727c88bea45049464947bad..e2b085f76a8e951886d668915306667a7de1b5f2 100644 (file)
@@ -22,31 +22,7 @@ UnifracUnweightedCommand::UnifracUnweightedCommand() {
                openOutputFile(sumFile, outSum);
                distFile = globaldata->getTreeFile() + ".uwdistrib";
                openOutputFile(distFile, outDist);
-
-               //if the user has not entered specific groups to analyze then do them all
-               if (globaldata->Groups.size() != 0) {
-                       //check that groups are valid
-                       for (int i = 0; i < globaldata->Groups.size(); i++) {
-                               if (tmap->isValidGroup(globaldata->Groups[i]) != true) {
-                                       cout << globaldata->Groups[i] << " is not a valid group, and will be disregarded." << endl;
-                                       // erase the invalid group from globaldata->Groups
-                                       globaldata->Groups.erase (globaldata->Groups.begin()+i);
-                               }
-                       }
-                       
-                       //if the user only entered invalid groups
-                       if (globaldata->Groups.size() == 0) { 
-                               cout << "When using the groups parameter you must have at least 1 valid group. I will run the command using all the groups in your groupfile." << endl; 
-                               for (int i = 0; i < tmap->namesOfGroups.size(); i++) {
-                                       globaldata->Groups.push_back(tmap->namesOfGroups[i]);
-                               }
-                       }
-               }else {
-                       for (int i = 0; i < tmap->namesOfGroups.size(); i++) {
-                               globaldata->Groups.push_back(tmap->namesOfGroups[i]);
-                       }
-               }
-               
+               setGroups(); //sets users groups to analyze
                convert(globaldata->getIters(), iters);  //how many random trees to generate
                unweighted = new Unweighted(tmap);
 
@@ -70,11 +46,18 @@ int UnifracUnweightedCommand::execute() {
                
                //format output
                outDist.setf(ios::fixed, ios::floatfield); outDist.setf(ios::showpoint);
+               
+               outDist << "Groups Used ";
+               for (int m = 0; m < globaldata->Groups.size(); m++) {
+                       outDist << globaldata->Groups[m] << " ";
+               }
+               outDist << endl;
+               
                outDist << "Tree#" << '\t' << "Iter" << '\t' << "UWScore" << endl;
                
                //create new tree with same num nodes and leaves as users
                randT = new Tree();
-                       
+                               
                //get pscores for users trees
                for (int i = 0; i < T.size(); i++) {
                        cout << "Processing tree " << i+1 << endl;
@@ -183,6 +166,12 @@ void UnifracUnweightedCommand::printUnweightedFile() {
        try {
                //column headers
                
+               out << "Groups Used ";
+               for (int m = 0; m < globaldata->Groups.size(); m++) {
+                       out << globaldata->Groups[m] << " ";
+               }
+               out << endl;
+
                out << "Score" << '\t' << "UserFreq" << '\t' << "UserCumul" << '\t' << "RandFreq" << '\t' << "RandCumul" << endl;
                                
                //format output
@@ -210,6 +199,13 @@ void UnifracUnweightedCommand::printUnweightedFile() {
 void UnifracUnweightedCommand::printUWSummaryFile() {
        try {
                //column headers
+               
+               outSum << "Groups Used ";
+               for (int m = 0; m < globaldata->Groups.size(); m++) {
+                       outSum << globaldata->Groups[m] << " ";
+               }
+               outSum << endl;
+
                outSum << "Tree#" << '\t'  <<  "UWScore" << '\t' << '\t' << "UWSig" <<  endl;
                
                //format output
@@ -255,4 +251,44 @@ void UnifracUnweightedCommand::saveRandomScores() {
        }
 }
 
-/***********************************************************/
\ No newline at end of file
+/***********************************************************/
+
+void UnifracUnweightedCommand::setGroups() {
+       try {
+               //if the user has not entered specific groups to analyze then do them all
+               if (globaldata->Groups.size() != 0) {
+                       //check that groups are valid
+                       for (int i = 0; i < globaldata->Groups.size(); i++) {
+                               if (tmap->isValidGroup(globaldata->Groups[i]) != true) {
+                                       cout << globaldata->Groups[i] << " is not a valid group, and will be disregarded." << endl;
+                                       // erase the invalid group from globaldata->Groups
+                                       globaldata->Groups.erase (globaldata->Groups.begin()+i);
+                               }
+                       }
+                       
+                       //if the user only entered invalid groups
+                       if (globaldata->Groups.size() == 0) { 
+                               cout << "When using the groups parameter you must have at least 1 valid group. I will run the command using all the groups in your groupfile." << endl; 
+                               for (int i = 0; i < tmap->namesOfGroups.size(); i++) {
+                                       globaldata->Groups.push_back(tmap->namesOfGroups[i]);
+                               }
+                       }
+                                       
+               }else {
+                       for (int i = 0; i < tmap->namesOfGroups.size(); i++) {
+                               globaldata->Groups.push_back(tmap->namesOfGroups[i]);
+                       }
+               }
+       }
+       catch(exception& e) {
+               cout << "Standard Error: " << e.what() << " has occurred in the UnifracUnweightedCommand class Function setGroups. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
+               exit(1);
+       }
+       catch(...) {
+               cout << "An unknown error has occurred in the UnifracUnweightedCommand class function setGroups. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
+               exit(1);
+       }               
+
+}
+/*****************************************************************/
+
index 9fe164298b0c9191674a01b4a4fd72fc4fed37ae..7bf53386a3100f0c881a1f937ad18b887132ca00 100644 (file)
@@ -50,7 +50,8 @@ class UnifracUnweightedCommand : public Command {
                
                void printUWSummaryFile();
                void printUnweightedFile();
-               void saveRandomScores();    
+               void saveRandomScores(); 
+               void setGroups();   
                
 };
 
index 4416e79afd680059b8958cab6c28e035b3633d1c..1a2573b6d5738055ddd7d69448a321369d72bb53 100644 (file)
@@ -17,56 +17,81 @@ EstOutput Unweighted::getValues(Tree* t) {
                
                //clear out old values
                data.resize(1,0); 
-               penalty.resize(t->getNumLeaves(), 0);
                
-               map<string,double> unique;  //group, total of all branch lengths of nodes with that group.
-               double shared = 0.0000;
-               double UW=0.0000;
+               float UniqueBL=0.0000;  //a branch length is unique if it's chidren are from the same group
+               float totalBL = 0.00;   //all branch lengths
+               float UW = 0.00;                //Unweighted Value = UniqueBL / totalBL;
                
-               //add up the branch lengths for each group. 
-               for(int i=0;i<t->getNumLeaves();i++){
-                       if(t->tree[i].pGroups.size() > 0){
-                               unique[t->tree[i].pGroups.begin()->first] += t->tree[i].getBranchLength();
-                       }
-               }
-               
-               //for each non-leaf node
+               map<string, int>::iterator it;  //iterator to traverse pgroups
+               map<string, int> copyLCpcount;
+               map<string, int> copyRCpcount;
+               map<string, int> copyIpcount;
+       
                for(int i=t->getNumLeaves();i<t->getNumNodes();i++){
                
                        int lc = t->tree[i].getLChild();  //lc = vector index of left child
                        int rc = t->tree[i].getRChild();  //rc = vector index of right child
                        
-                       //get penalty values
-                       if(t->tree[rc].pGroups.size() == 0 || t->tree[lc].pGroups.size() == 0){
-                               penalty.push_back(penalty[t->tree[rc].getIndex()]+penalty[t->tree[lc].getIndex()]);
-                       }
-                       else if(t->tree[i].pGroups.size() > t->tree[rc].pGroups.size() || t->tree[i].pGroups.size() > t->tree[lc].pGroups.size()){
-                               penalty.push_back(penalty[t->tree[rc].getIndex()]+penalty[t->tree[lc].getIndex()]+1);
+                       /**********************************************************************/
+                       //This section adds in all lengths that are non leaf
+                       
+                       //copy left childs pGroups and remove groups that the user doesn't want
+                       copyIpcount = t->tree[i].pcount;
+                       for (it = copyIpcount.begin(); it != copyIpcount.end(); it++) {
+                               if (inUsersGroups(it->first, globaldata->Groups) != true) {     copyIpcount.erase(it->first);   }
                        }
-                       else{
-                               penalty.push_back(penalty[t->tree[rc].getIndex()]+penalty[t->tree[lc].getIndex()]);
+
+                       //copy left childs pGroups and remove groups that the user doesn't want
+                       copyLCpcount = t->tree[lc].pcount;
+                       for (it = copyLCpcount.begin(); it != copyLCpcount.end(); it++) {
+                               if (inUsersGroups(it->first, globaldata->Groups) != true) {     copyLCpcount.erase(it->first);  }
                        }
 
-                       //not sure when this would ever be true??? if your parent is root could be, but pGroups.size() should never be 0.
-                       if(t->tree[i].getParent() == -1 && (t->tree[lc].pGroups.size() == 0 || t->tree[rc].pGroups.size() == 0)){
-                               shared -= 1; 
+                       //copy right childs pGroups and remove groups that the user doesn't want
+                       copyRCpcount = t->tree[rc].pcount;
+                       for (it = copyRCpcount.begin(); it != copyRCpcount.end(); it++) {
+                               if (inUsersGroups(it->first, globaldata->Groups) != true) {     copyRCpcount.erase(it->first);  }
                        }
-                       else if(penalty[i] != 0 && t->tree[i].pGroups.size() != 0){
-                               shared += t->tree[i].getBranchLength();
+       
+                       //if i's children are from the same group and i has a BL then add i's length to unique
+                       //if copyRCpcount.size() = 0 && copyLCpcount.size() = 0 they are from a branch that is entirely from a group the user doesn't want
+                       if ((copyRCpcount.size() == 0) && (copyLCpcount.size() == 0)) { }
+                       else {
+                               if ((copyRCpcount == copyLCpcount) && (t->tree[i].getBranchLength() != -1)) {  UniqueBL += t->tree[i].getBranchLength();        }
+                               //if either childs groups = 0 then all of there groups were not valid making the parent unique
+                               else if (((copyRCpcount.size() == 0) || (copyLCpcount.size() == 0)) && (t->tree[i].getBranchLength() != -1)) {  UniqueBL += t->tree[i].getBranchLength();       }
                        }
-                       else if( t->tree[i].pGroups.size() != 0){
-                               unique[t->tree[i].pGroups.begin()->first] += t->tree[i].getBranchLength();          
+                       
+                       //add i's BL to total if it is from the groups the user wants
+                       if ((t->tree[i].getBranchLength() != -1) && (copyIpcount.size() != 0)) {  
+                               totalBL += t->tree[i].getBranchLength(); 
                        }
-               }
-    
-               map<string,double>::iterator pos;
-               for(pos=unique.begin();pos!=unique.end();pos++){
-                       if((pos->first!="xxx") && (inUsersGroups(pos->first))){     
-                               UW += unique[pos->first];
+                       
+                       /**********************************************************************/
+                       //This section adds in all lengths that are leaf
+                       
+                       //if i's chidren are leaves
+                       if (t->tree[rc].getRChild() == -1) {
+                               //if rc is a valid group and rc has a BL
+                               if ((inUsersGroups(t->tree[rc].getGroup(), globaldata->Groups) == true) && (t->tree[rc].getBranchLength() != -1)) {
+                                       UniqueBL += t->tree[rc].getBranchLength();
+                                       totalBL += t->tree[rc].getBranchLength(); 
+                               }
+                       }
+                       
+                       if (t->tree[lc].getLChild() == -1) {
+                               //if lc is a valid group and lc has a BL
+                               if ((inUsersGroups(t->tree[lc].getGroup(), globaldata->Groups) == true) && (t->tree[lc].getBranchLength() != -1)) {
+                                       UniqueBL += t->tree[lc].getBranchLength();
+                                       totalBL += t->tree[lc].getBranchLength(); 
+                               }
                        }
+                       
+                       /**********************************************************************/
+               
                }
-       
-               UW /= (UW + shared);
+               
+               UW = (UniqueBL / totalBL);  
        
                if (isnan(UW) || isinf(UW)) { UW = 0; }
        
@@ -86,20 +111,3 @@ EstOutput Unweighted::getValues(Tree* t) {
 
 }
 
-/**************************************************************************************************/
-bool Unweighted::inUsersGroups(string groupname) {
-       try {
-               for (int i = 0; i < globaldata->Groups.size(); i++) {
-                       if (groupname == globaldata->Groups[i]) { return true; }
-               }
-               return false;
-       }
-       catch(exception& e) {
-               cout << "Standard Error: " << e.what() << " has occurred in the Unweighted class Function inUsersGroups. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
-               exit(1);
-       }
-       catch(...) {
-               cout << "An unknown error has occurred in the Unweighted class function inUsersGroups. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
-               exit(1);
-       }
-}
\ No newline at end of file
index e3bf3a8fc240648e4a38f48ced539845d3648706..b6f225f07392fde58d7afd46536562ee5514b104 100644 (file)
@@ -26,9 +26,7 @@ class Unweighted : public TreeCalculator  {
        private:
                GlobalData* globaldata;
                EstOutput data;
-               vector<int> penalty;
                TreeMap* tmap;
-               bool inUsersGroups(string);
 
 };
 
index cc6c4900d459828e5fbd2a9336754e589e160266..90d10c1cc4e96f34561c37f81738015158957488 100644 (file)
@@ -4,6 +4,7 @@
 using namespace std;
 
 #include "mothur.h"
+#include "treemap.h"
 
 typedef unsigned long long ull;
 
@@ -333,8 +334,23 @@ inline void splitAtEquals(string& key, string& value){
        }
 
 }
-/*******************************************************/
-
+/**************************************************************************************************/
 
+inline bool inUsersGroups(string groupname, vector<string> Groups) {
+       try {
+               for (int i = 0; i < Groups.size(); i++) {
+                       if (groupname == Groups[i]) { return true; }
+               }
+               return false;
+       }
+       catch(exception& e) {
+               cout << "Standard Error: " << e.what() << " has occurred in the utilities class Function inUsersGroups. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
+               exit(1);
+       }
+       catch(...) {
+               cout << "An unknown error has occurred in the utilities class function inUsersGroups. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
+               exit(1);
+       }
+}
 
 #endif
index fca5c8e07fc50d9be08842e0787854431a6df4a9..590946e74101ff7cdc37779e86504b75567937c5 100644 (file)
@@ -34,12 +34,6 @@ ValidParameters::ValidParameters() {
                parameters["iters"]                             = "iters"; 
                parameters["jumble"]                    = "jumble"; 
                parameters["freq"]                              = "freq"; 
-               parameters["single"]                    = "single"; 
-               parameters["rarefaction"]               = "rarefaction"; 
-               parameters["sharedrarefaction"] = "sharedrarefaction";
-               parameters["shared"]                    = "shared"; 
-               parameters["summary"]                   = "summary"; 
-               parameters["sharedsummary"]             = "sharedsummary";
                parameters["abund"]             = "abund";
                parameters["random"]                    = "random";
                parameters["groups"]                    = "groups";