From d97b619c4297b1274c754d73a64792ba656b0a79 Mon Sep 17 00:00:00 2001 From: westcott Date: Fri, 20 Feb 2009 17:34:55 +0000 Subject: [PATCH] fixed unweighted calculator --- errorchecking.cpp | 45 +++++++------- errorchecking.h | 10 ++-- parsimony.cpp | 24 +------- parsimony.h | 1 - parsimonycommand.cpp | 67 +++++++++++++-------- parsimonycommand.h | 1 + tree.cpp | 77 ++++++++++++++++++++---- tree.h | 1 + unifracunweightedcommand.cpp | 90 +++++++++++++++++++--------- unifracunweightedcommand.h | 3 +- unweighted.cpp | 112 +++++++++++++++++++---------------- unweighted.h | 2 - utilities.hpp | 20 ++++++- validparameter.cpp | 6 -- 14 files changed, 283 insertions(+), 176 deletions(-) diff --git a/errorchecking.cpp b/errorchecking.cpp index dccb679..ff94c8f 100644 --- a/errorchecking.cpp +++ b/errorchecking.cpp @@ -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 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 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) { diff --git a/errorchecking.h b/errorchecking.h index 8bf6857..50cc005 100644 --- a/errorchecking.h +++ b/errorchecking.h @@ -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 sharedGroups; - map > commandParameters; - map > intParams; - double piSent; + // map > commandParameters; + /// map > intParams; + // double piSent; }; #endif diff --git a/parsimony.cpp b/parsimony.cpp index 2a673b0..6239f75 100644 --- a/parsimony.cpp +++ b/parsimony.cpp @@ -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); - } -} -/**************************************************************************************************/ diff --git a/parsimony.h b/parsimony.h index 8955ef7..ebe3d8d 100644 --- a/parsimony.h +++ b/parsimony.h @@ -28,7 +28,6 @@ class Parsimony : public TreeCalculator { GlobalData* globaldata; EstOutput data; TreeMap* tmap; - bool inUsersGroups(string); map::iterator it; }; diff --git a/parsimonycommand.cpp b/parsimonycommand.cpp index fe5a1dd..554d24b 100644 --- a/parsimonycommand.cpp +++ b/parsimonycommand.cpp @@ -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); + } + +} +/*****************************************************************/ + + diff --git a/parsimonycommand.h b/parsimonycommand.h index bd52fa6..326c0c9 100644 --- a/parsimonycommand.h +++ b/parsimonycommand.h @@ -51,6 +51,7 @@ class ParsimonyCommand : public Command { void printParsimonyFile(); void printUSummaryFile(); void getUserInput(); + void setGroups(); }; diff --git a/tree.cpp b/tree.cpp index 94ab326..2ad46f0 100644 --- a/tree.cpp +++ b/tree.cpp @@ -231,33 +231,51 @@ map 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 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 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 gcount_hold = tree[z].pcount; - tree[z].pcount = (tree[i].pcount); - tree[i].pcount = (gcount_hold); + map 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 d7dcc3d..7157a6f 100644 --- 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 diff --git a/unifracunweightedcommand.cpp b/unifracunweightedcommand.cpp index bf2129c..e2b085f 100644 --- a/unifracunweightedcommand.cpp +++ b/unifracunweightedcommand.cpp @@ -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); + } + +} +/*****************************************************************/ + diff --git a/unifracunweightedcommand.h b/unifracunweightedcommand.h index 9fe1642..7bf5338 100644 --- a/unifracunweightedcommand.h +++ b/unifracunweightedcommand.h @@ -50,7 +50,8 @@ class UnifracUnweightedCommand : public Command { void printUWSummaryFile(); void printUnweightedFile(); - void saveRandomScores(); + void saveRandomScores(); + void setGroups(); }; diff --git a/unweighted.cpp b/unweighted.cpp index 4416e79..1a2573b 100644 --- a/unweighted.cpp +++ b/unweighted.cpp @@ -17,56 +17,81 @@ EstOutput Unweighted::getValues(Tree* t) { //clear out old values data.resize(1,0); - penalty.resize(t->getNumLeaves(), 0); - map 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;igetNumLeaves();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::iterator it; //iterator to traverse pgroups + map copyLCpcount; + map copyRCpcount; + map copyIpcount; + for(int i=t->getNumLeaves();igetNumNodes();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::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 diff --git a/unweighted.h b/unweighted.h index e3bf3a8..b6f225f 100644 --- a/unweighted.h +++ b/unweighted.h @@ -26,9 +26,7 @@ class Unweighted : public TreeCalculator { private: GlobalData* globaldata; EstOutput data; - vector penalty; TreeMap* tmap; - bool inUsersGroups(string); }; diff --git a/utilities.hpp b/utilities.hpp index cc6c490..90d10c1 100644 --- a/utilities.hpp +++ b/utilities.hpp @@ -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 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 diff --git a/validparameter.cpp b/validparameter.cpp index fca5c8e..590946e 100644 --- a/validparameter.cpp +++ b/validparameter.cpp @@ -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"; -- 2.39.2