2 * phylodiversitycommand.cpp
5 * Created by westcott on 4/30/10.
6 * Copyright 2010 Schloss Lab. All rights reserved.
10 #include "phylodiversitycommand.h"
11 #include "treereader.h"
13 //**********************************************************************************************************************
14 vector<string> PhyloDiversityCommand::setParameters(){
17 CommandParameter ptree("tree", "InputTypes", "", "", "none", "none", "none",false,true); parameters.push_back(ptree);
18 CommandParameter pgroup("group", "InputTypes", "", "", "none", "none", "none",false,true); parameters.push_back(pgroup);
19 CommandParameter pname("name", "InputTypes", "", "", "none", "none", "none",false,false); parameters.push_back(pname);
20 CommandParameter pgroups("groups", "String", "", "", "", "", "",false,false); parameters.push_back(pgroups);
21 CommandParameter piters("iters", "Number", "", "1000", "", "", "",false,false); parameters.push_back(piters);
22 CommandParameter pfreq("freq", "Number", "", "100", "", "", "",false,false); parameters.push_back(pfreq);
23 CommandParameter pprocessors("processors", "Number", "", "1", "", "", "",false,false); parameters.push_back(pprocessors);
24 CommandParameter prarefy("rarefy", "Boolean", "", "F", "", "", "",false,false); parameters.push_back(prarefy);
25 CommandParameter psummary("summary", "Boolean", "", "T", "", "", "",false,false); parameters.push_back(psummary);
26 CommandParameter pcollect("collect", "Boolean", "", "F", "", "", "",false,false); parameters.push_back(pcollect);
27 CommandParameter pscale("scale", "Boolean", "", "F", "", "", "",false,false); parameters.push_back(pscale);
28 CommandParameter pinputdir("inputdir", "String", "", "", "", "", "",false,false); parameters.push_back(pinputdir);
29 CommandParameter poutputdir("outputdir", "String", "", "", "", "", "",false,false); parameters.push_back(poutputdir);
31 vector<string> myArray;
32 for (int i = 0; i < parameters.size(); i++) { myArray.push_back(parameters[i].name); }
36 m->errorOut(e, "PhyloDiversityCommand", "setParameters");
40 //**********************************************************************************************************************
41 string PhyloDiversityCommand::getHelpString(){
43 string helpString = "";
44 helpString += "The phylo.diversity command parameters are tree, group, name, groups, iters, freq, processors, scale, rarefy, collect and summary. tree and group are required, unless you have valid current files.\n";
45 helpString += "The groups parameter allows you to specify which of the groups in your groupfile you would like analyzed. The group names are separated by dashes. By default all groups are used.\n";
46 helpString += "The iters parameter allows you to specify the number of randomizations to preform, by default iters=1000, if you set rarefy to true.\n";
47 helpString += "The freq parameter is used indicate when to output your data, by default it is set to 100. But you can set it to a percentage of the number of sequence. For example freq=0.10, means 10%. \n";
48 helpString += "The scale parameter is used indicate that you want your output scaled to the number of sequences sampled, default = false. \n";
49 helpString += "The rarefy parameter allows you to create a rarefaction curve. The default is false.\n";
50 helpString += "The collect parameter allows you to create a collectors curve. The default is false.\n";
51 helpString += "The summary parameter allows you to create a .summary file. The default is true.\n";
52 helpString += "The processors parameter allows you to specify the number of processors to use. The default is 1.\n";
53 helpString += "The phylo.diversity command should be in the following format: phylo.diversity(groups=yourGroups, rarefy=yourRarefy, iters=yourIters).\n";
54 helpString += "Example phylo.diversity(groups=A-B-C, rarefy=T, iters=500).\n";
55 helpString += "The phylo.diversity command output two files: .phylo.diversity and if rarefy=T, .rarefaction.\n";
56 helpString += "Note: No spaces between parameter labels (i.e. groups), '=' and parameters (i.e.yourGroups).\n";
60 m->errorOut(e, "PhyloDiversityCommand", "getHelpString");
64 //**********************************************************************************************************************
65 string PhyloDiversityCommand::getOutputFileNameTag(string type, string inputName=""){
67 string outputFileName = "";
68 map<string, vector<string> >::iterator it;
70 //is this a type this command creates
71 it = outputTypes.find(type);
72 if (it == outputTypes.end()) { m->mothurOut("[ERROR]: this command doesn't create a " + type + " output file.\n"); }
74 if (type == "phylodiv") { outputFileName = "phylodiv"; }
75 else if (type == "rarefy") { outputFileName = "phylodiv.rarefaction"; }
76 else if (type == "summary") { outputFileName = "phylodiv.summary"; }
77 else { m->mothurOut("[ERROR]: No definition for type " + type + " output file tag.\n"); m->control_pressed = true; }
79 return outputFileName;
82 m->errorOut(e, "PhyloDiversityCommand", "getOutputFileNameTag");
87 //**********************************************************************************************************************
88 PhyloDiversityCommand::PhyloDiversityCommand(){
90 abort = true; calledHelp = true;
92 vector<string> tempOutNames;
93 outputTypes["phylodiv"] = tempOutNames;
94 outputTypes["rarefy"] = tempOutNames;
95 outputTypes["summary"] = tempOutNames;
98 m->errorOut(e, "PhyloDiversityCommand", "PhyloDiversityCommand");
102 //**********************************************************************************************************************
103 PhyloDiversityCommand::PhyloDiversityCommand(string option) {
105 abort = false; calledHelp = false;
107 //allow user to run help
108 if(option == "help") { help(); abort = true; calledHelp = true; }
109 else if(option == "citation") { citation(); abort = true; calledHelp = true;}
112 vector<string> myArray = setParameters();;
114 OptionParser parser(option);
115 map<string,string> parameters = parser.getParameters();
116 map<string,string>::iterator it;
118 ValidParameters validParameter;
120 //check to make sure all parameters are valid for command
121 for (map<string,string>::iterator it = parameters.begin(); it != parameters.end(); it++) {
122 if (validParameter.isValidParameter(it->first, myArray, it->second) != true) { abort = true; }
125 //initialize outputTypes
126 vector<string> tempOutNames;
127 outputTypes["phylodiv"] = tempOutNames;
128 outputTypes["rarefy"] = tempOutNames;
129 outputTypes["summary"] = tempOutNames;
131 //if the user changes the input directory command factory will send this info to us in the output parameter
132 string inputDir = validParameter.validFile(parameters, "inputdir", false);
133 if (inputDir == "not found"){ inputDir = ""; }
136 it = parameters.find("tree");
137 //user has given a template file
138 if(it != parameters.end()){
139 path = m->hasPath(it->second);
140 //if the user has not given a path then, add inputdir. else leave path alone.
141 if (path == "") { parameters["tree"] = inputDir + it->second; }
144 it = parameters.find("group");
145 //user has given a template file
146 if(it != parameters.end()){
147 path = m->hasPath(it->second);
148 //if the user has not given a path then, add inputdir. else leave path alone.
149 if (path == "") { parameters["group"] = inputDir + it->second; }
152 it = parameters.find("name");
153 //user has given a template file
154 if(it != parameters.end()){
155 path = m->hasPath(it->second);
156 //if the user has not given a path then, add inputdir. else leave path alone.
157 if (path == "") { parameters["name"] = inputDir + it->second; }
161 //check for required parameters
162 treefile = validParameter.validFile(parameters, "tree", true);
163 if (treefile == "not open") { treefile = ""; abort = true; }
164 else if (treefile == "not found") {
165 //if there is a current design file, use it
166 treefile = m->getTreeFile();
167 if (treefile != "") { m->mothurOut("Using " + treefile + " as input file for the tree parameter."); m->mothurOutEndLine(); }
168 else { m->mothurOut("You have no current tree file and the tree parameter is required."); m->mothurOutEndLine(); abort = true; }
169 }else { m->setTreeFile(treefile); }
171 //check for required parameters
172 groupfile = validParameter.validFile(parameters, "group", true);
173 if (groupfile == "not open") { groupfile = ""; abort = true; }
174 else if (groupfile == "not found") { groupfile = ""; }
175 else { m->setGroupFile(groupfile); }
177 namefile = validParameter.validFile(parameters, "name", true);
178 if (namefile == "not open") { namefile = ""; abort = true; }
179 else if (namefile == "not found") { namefile = ""; }
180 else { m->setNameFile(namefile); }
182 outputDir = validParameter.validFile(parameters, "outputdir", false); if (outputDir == "not found"){ outputDir = m->hasPath(treefile); }
185 temp = validParameter.validFile(parameters, "freq", false); if (temp == "not found") { temp = "100"; }
186 m->mothurConvert(temp, freq);
188 temp = validParameter.validFile(parameters, "iters", false); if (temp == "not found") { temp = "1000"; }
189 m->mothurConvert(temp, iters);
191 temp = validParameter.validFile(parameters, "rarefy", false); if (temp == "not found") { temp = "F"; }
192 rarefy = m->isTrue(temp);
193 if (!rarefy) { iters = 1; }
195 temp = validParameter.validFile(parameters, "summary", false); if (temp == "not found") { temp = "T"; }
196 summary = m->isTrue(temp);
198 temp = validParameter.validFile(parameters, "scale", false); if (temp == "not found") { temp = "F"; }
199 scale = m->isTrue(temp);
201 temp = validParameter.validFile(parameters, "collect", false); if (temp == "not found") { temp = "F"; }
202 collect = m->isTrue(temp);
204 temp = validParameter.validFile(parameters, "processors", false); if (temp == "not found"){ temp = m->getProcessors(); }
205 m->setProcessors(temp);
206 m->mothurConvert(temp, processors);
208 groups = validParameter.validFile(parameters, "groups", false);
209 if (groups == "not found") { groups = ""; }
211 m->splitAtDash(groups, Groups);
212 m->setGroups(Groups);
215 if ((!collect) && (!rarefy) && (!summary)) { m->mothurOut("No outputs selected. You must set either collect, rarefy or summary to true, summary=T by default."); m->mothurOutEndLine(); abort=true; }
217 if (namefile == "") {
218 vector<string> files; files.push_back(treefile);
219 parser.getNameFile(files);
224 catch(exception& e) {
225 m->errorOut(e, "PhyloDiversityCommand", "PhyloDiversityCommand");
229 //**********************************************************************************************************************
231 int PhyloDiversityCommand::execute(){
234 if (abort == true) { if (calledHelp) { return 0; } return 2; }
236 int start = time(NULL);
238 m->setTreeFile(treefile);
239 TreeReader* reader = new TreeReader(treefile, groupfile, namefile);
240 vector<Tree*> trees = reader->getTrees();
241 tmap = trees[0]->getTreeMap();
245 vector<string> mGroups = m->getGroups();
246 vector<string> tGroups = tmap->getNamesOfGroups();
247 util.setGroups(mGroups, tGroups, "phylo.diversity"); //sets the groups the user wants to analyze
249 //incase the user had some mismatches between the tree and group files we don't want group xxx to be analyzed
250 for (int i = 0; i < mGroups.size(); i++) { if (mGroups[i] == "xxx") { mGroups.erase(mGroups.begin()+i); break; } }
251 m->setGroups(mGroups);
253 vector<string> outputNames;
255 //for each of the users trees
256 for(int i = 0; i < trees.size(); i++) {
258 if (m->control_pressed) { delete tmap; for (int j = 0; j < trees.size(); j++) { delete trees[j]; } for (int j = 0; j < outputNames.size(); j++) { m->mothurRemove(outputNames[j]); } return 0; }
260 ofstream outSum, outRare, outCollect;
261 string outSumFile = outputDir + m->getRootName(m->getSimpleName(treefile)) + toString(i+1) + "." + getOutputFileNameTag("summary");
262 string outRareFile = outputDir + m->getRootName(m->getSimpleName(treefile)) + toString(i+1) + "." + getOutputFileNameTag("rarefy");
263 string outCollectFile = outputDir + m->getRootName(m->getSimpleName(treefile)) + toString(i+1) + "." + getOutputFileNameTag("phylodiv");
265 if (summary) { m->openOutputFile(outSumFile, outSum); outputNames.push_back(outSumFile); outputTypes["summary"].push_back(outSumFile); }
266 if (rarefy) { m->openOutputFile(outRareFile, outRare); outputNames.push_back(outRareFile); outputTypes["rarefy"].push_back(outRareFile); }
267 if (collect) { m->openOutputFile(outCollectFile, outCollect); outputNames.push_back(outCollectFile); outputTypes["phylodiv"].push_back(outCollectFile); }
269 int numLeafNodes = trees[i]->getNumLeaves();
271 //create a vector containing indexes of leaf nodes, randomize it, select nodes to send to calculator
272 vector<int> randomLeaf;
273 for (int j = 0; j < numLeafNodes; j++) {
274 if (m->inUsersGroups(trees[i]->tree[j].getGroup(), mGroups) == true) { //is this a node from the group the user selected.
275 randomLeaf.push_back(j);
279 numLeafNodes = randomLeaf.size(); //reset the number of leaf nodes you are using
281 //each group, each sampling, if no rarefy iters = 1;
282 map<string, vector<float> > diversity;
284 //each group, each sampling, if no rarefy iters = 1;
285 map<string, vector<float> > sumDiversity;
287 //find largest group total
288 int largestGroup = 0;
289 for (int j = 0; j < mGroups.size(); j++) {
290 if (tmap->seqsPerGroup[mGroups[j]] > largestGroup) { largestGroup = tmap->seqsPerGroup[mGroups[j]]; }
292 //initialize diversity
293 diversity[mGroups[j]].resize(tmap->seqsPerGroup[mGroups[j]]+1, 0.0); //numSampled
296 //initialize sumDiversity
297 sumDiversity[mGroups[j]].resize(tmap->seqsPerGroup[mGroups[j]]+1, 0.0);
300 //convert freq percentage to number
302 if (freq < 1.0) { increment = largestGroup * freq;
303 }else { increment = freq; }
305 //initialize sampling spots
306 set<int> numSampledList;
307 for(int k = 1; k <= largestGroup; k++){ if((k == 1) || (k % increment == 0)){ numSampledList.insert(k); } }
308 if(largestGroup % increment != 0){ numSampledList.insert(largestGroup); }
310 //add other groups ending points
311 for (int j = 0; j < mGroups.size(); j++) {
312 if (numSampledList.count(diversity[mGroups[j]].size()-1) == 0) { numSampledList.insert(diversity[mGroups[j]].size()-1); }
315 #if defined (__APPLE__) || (__MACH__) || (linux) || (__linux) || (__linux__) || (__unix__) || (__unix)
317 driver(trees[i], diversity, sumDiversity, iters, increment, randomLeaf, numSampledList, outCollect, outSum, true);
320 vector<int> procIters;
322 int numItersPerProcessor = iters / processors;
324 //divide iters between processes
325 for (int h = 0; h < processors; h++) {
326 if(h == processors - 1){
327 numItersPerProcessor = iters - h * numItersPerProcessor;
329 procIters.push_back(numItersPerProcessor);
332 createProcesses(procIters, trees[i], diversity, sumDiversity, iters, increment, randomLeaf, numSampledList, outCollect, outSum);
334 }else{ //no need to paralellize if you dont want to rarefy
335 driver(trees[i], diversity, sumDiversity, iters, increment, randomLeaf, numSampledList, outCollect, outSum, true);
340 driver(trees[i], diversity, sumDiversity, iters, increment, randomLeaf, numSampledList, outCollect, outSum, true);
343 if (rarefy) { printData(numSampledList, sumDiversity, outRare, iters); }
347 if (m->control_pressed) { for (int i = 0; i < outputNames.size(); i++) { m->mothurRemove(outputNames[i]); } return 0; }
349 m->mothurOut("It took " + toString(time(NULL) - start) + " secs to run phylo.diversity."); m->mothurOutEndLine();
352 m->mothurOutEndLine();
353 m->mothurOut("Output File Names: "); m->mothurOutEndLine();
354 for (int i = 0; i < outputNames.size(); i++) { m->mothurOut(outputNames[i]); m->mothurOutEndLine(); }
355 m->mothurOutEndLine();
360 catch(exception& e) {
361 m->errorOut(e, "PhyloDiversityCommand", "execute");
365 //**********************************************************************************************************************
366 int PhyloDiversityCommand::createProcesses(vector<int>& procIters, Tree* t, map< string, vector<float> >& div, map<string, vector<float> >& sumDiv, int numIters, int increment, vector<int>& randomLeaf, set<int>& numSampledList, ofstream& outCollect, ofstream& outSum){
368 #if defined (__APPLE__) || (__MACH__) || (linux) || (__linux) || (__linux__) || (__unix__) || (__unix)
371 vector<int> processIDS;
372 map< string, vector<float> >::iterator itSum;
374 //loop through and create all the processes you want
375 while (process != processors) {
379 processIDS.push_back(pid); //create map from line number to pid so you can append files in correct order later
382 driver(t, div, sumDiv, procIters[process], increment, randomLeaf, numSampledList, outCollect, outSum, false);
384 string outTemp = outputDir + toString(getpid()) + ".sumDiv.temp";
386 m->openOutputFile(outTemp, out);
388 //output the sumDIversity
389 for (itSum = sumDiv.begin(); itSum != sumDiv.end(); itSum++) {
390 out << itSum->first << '\t' << (itSum->second).size() << '\t';
391 for (int k = 0; k < (itSum->second).size(); k++) {
392 out << (itSum->second)[k] << '\t';
401 m->mothurOut("[ERROR]: unable to spawn the necessary processes."); m->mothurOutEndLine();
402 for (int i = 0; i < processIDS.size(); i++) { kill (processIDS[i], SIGINT); }
407 driver(t, div, sumDiv, procIters[0], increment, randomLeaf, numSampledList, outCollect, outSum, true);
409 //force parent to wait until all the processes are done
410 for (int i=0;i<(processors-1);i++) {
411 int temp = processIDS[i];
415 //get data created by processes
416 for (int i=0;i<(processors-1);i++) {
418 //input the sumDIversity
419 string inTemp = outputDir + toString(processIDS[i]) + ".sumDiv.temp";
421 m->openInputFile(inTemp, in);
423 //output the sumDIversity
424 for (int j = 0; j < sumDiv.size(); j++) {
428 in >> group >> size; m->gobble(in);
430 for (int k = 0; k < size; k++) {
434 sumDiv[group][k] += tempVal;
440 m->mothurRemove(inTemp);
448 catch(exception& e) {
449 m->errorOut(e, "PhyloDiversityCommand", "createProcesses");
453 //**********************************************************************************************************************
454 int PhyloDiversityCommand::driver(Tree* t, map< string, vector<float> >& div, map<string, vector<float> >& sumDiv, int numIters, int increment, vector<int>& randomLeaf, set<int>& numSampledList, ofstream& outCollect, ofstream& outSum, bool doSumCollect){
456 int numLeafNodes = randomLeaf.size();
457 vector<string> mGroups = m->getGroups();
459 map<string, int> rootForGroup = getRootForGroups(t); //maps groupName to root node in tree. "root" for group may not be the trees root and we don't want to include the extra branches.
461 for (int l = 0; l < numIters; l++) {
462 random_shuffle(randomLeaf.begin(), randomLeaf.end());
465 map<string, int> counts;
466 vector< map<string, bool> > countedBranch;
467 for (int i = 0; i < t->getNumNodes(); i++) {
468 map<string, bool> temp;
469 for (int j = 0; j < mGroups.size(); j++) { temp[mGroups[j]] = false; }
470 countedBranch.push_back(temp);
473 for (int j = 0; j < mGroups.size(); j++) { counts[mGroups[j]] = 0; }
475 for(int k = 0; k < numLeafNodes; k++){
477 if (m->control_pressed) { return 0; }
479 //calc branch length of randomLeaf k
480 vector<float> br = calcBranchLength(t, randomLeaf[k], countedBranch, rootForGroup);
482 //for each group in the groups update the total branch length accounting for the names file
483 vector<string> groups = t->tree[randomLeaf[k]].getGroup();
485 for (int j = 0; j < groups.size(); j++) {
487 if (m->inUsersGroups(groups[j], mGroups)) {
488 int numSeqsInGroupJ = 0;
489 map<string, int>::iterator it;
490 it = t->tree[randomLeaf[k]].pcount.find(groups[j]);
491 if (it != t->tree[randomLeaf[k]].pcount.end()) { //this leaf node contains seqs from group j
492 numSeqsInGroupJ = it->second;
495 if (numSeqsInGroupJ != 0) { div[groups[j]][(counts[groups[j]]+1)] = div[groups[j]][counts[groups[j]]] + br[j]; }
497 for (int s = (counts[groups[j]]+2); s <= (counts[groups[j]]+numSeqsInGroupJ); s++) {
498 div[groups[j]][s] = div[groups[j]][s-1]; //update counts, but don't add in redundant branch lengths
500 counts[groups[j]] += numSeqsInGroupJ;
506 //add this diversity to the sum
507 for (int j = 0; j < mGroups.size(); j++) {
508 for (int g = 0; g < div[mGroups[j]].size(); g++) {
509 sumDiv[mGroups[j]][g] += div[mGroups[j]][g];
514 if ((collect) && (l == 0) && doSumCollect) { printData(numSampledList, div, outCollect, 1); }
515 if ((summary) && (l == 0) && doSumCollect) { printSumData(div, outSum, 1); }
521 catch(exception& e) {
522 m->errorOut(e, "PhyloDiversityCommand", "driver");
527 //**********************************************************************************************************************
529 void PhyloDiversityCommand::printSumData(map< string, vector<float> >& div, ofstream& out, int numIters){
532 out << "Groups\tnumSampled\tphyloDiversity" << endl;
534 out.setf(ios::fixed, ios::floatfield); out.setf(ios::showpoint);
536 vector<string> mGroups = m->getGroups();
537 for (int j = 0; j < mGroups.size(); j++) {
538 int numSampled = (div[mGroups[j]].size()-1);
539 out << mGroups[j] << '\t' << numSampled << '\t';
543 if (scale) { score = (div[mGroups[j]][numSampled] / (float)numIters) / (float)numSampled; }
544 else { score = div[mGroups[j]][numSampled] / (float)numIters; }
546 out << setprecision(4) << score << endl;
552 catch(exception& e) {
553 m->errorOut(e, "PhyloDiversityCommand", "printSumData");
557 //**********************************************************************************************************************
559 void PhyloDiversityCommand::printData(set<int>& num, map< string, vector<float> >& div, ofstream& out, int numIters){
562 out << "numSampled\t";
563 vector<string> mGroups = m->getGroups();
564 for (int i = 0; i < mGroups.size(); i++) { out << mGroups[i] << '\t'; }
567 out.setf(ios::fixed, ios::floatfield); out.setf(ios::showpoint);
569 for (set<int>::iterator it = num.begin(); it != num.end(); it++) {
570 int numSampled = *it;
572 out << numSampled << '\t';
574 for (int j = 0; j < mGroups.size(); j++) {
575 if (numSampled < div[mGroups[j]].size()) {
577 if (scale) { score = (div[mGroups[j]][numSampled] / (float)numIters) / (float)numSampled; }
578 else { score = div[mGroups[j]][numSampled] / (float)numIters; }
580 out << setprecision(4) << score << '\t';
581 }else { out << "NA" << '\t'; }
589 catch(exception& e) {
590 m->errorOut(e, "PhyloDiversityCommand", "printData");
594 //**********************************************************************************************************************
595 //need a vector of floats one branch length for every group the node represents.
596 vector<float> PhyloDiversityCommand::calcBranchLength(Tree* t, int leaf, vector< map<string, bool> >& counted, map<string, int> roots){
599 //calc the branch length
600 //while you aren't at root
604 vector<string> groups = t->tree[leaf].getGroup();
605 sums.resize(groups.size(), 0.0);
609 if(t->tree[index].getBranchLength() != -1){
610 for (int k = 0; k < groups.size(); k++) {
611 sums[k] += abs(t->tree[index].getBranchLength());
616 index = t->tree[index].getParent();
618 //while you aren't at root
619 while(t->tree[index].getParent() != -1){
621 if (m->control_pressed) { return sums; }
623 for (int k = 0; k < groups.size(); k++) {
625 if (index >= roots[groups[k]]) { counted[index][groups[k]] = true; } //if you are at this groups "root", then say we are done
627 if (!counted[index][groups[k]]){ //if counted[index][groups[k] is true this groups has already added all br from here to root, so quit early
628 if (t->tree[index].getBranchLength() != -1) {
629 sums[k] += abs(t->tree[index].getBranchLength());
631 counted[index][groups[k]] = true;
634 index = t->tree[index].getParent();
640 catch(exception& e) {
641 m->errorOut(e, "PhyloDiversityCommand", "calcBranchLength");
645 //**********************************************************************************************************************
646 map<string, int> PhyloDiversityCommand::getRootForGroups(Tree* t){
648 map<string, int> roots; //maps group to root for group, may not be root of tree
649 map<string, bool> done;
651 //initialize root for all groups to -1
652 for (int k = 0; k < (t->getTreeMap())->getNamesOfGroups().size(); k++) { done[(t->getTreeMap())->getNamesOfGroups()[k]] = false; }
654 for (int i = 0; i < t->getNumLeaves(); i++) {
656 vector<string> groups = t->tree[i].getGroup();
658 int index = t->tree[i].getParent();
660 for (int j = 0; j < groups.size(); j++) {
662 if (done[groups[j]] == false) { //we haven't found the root for this group yet
664 done[groups[j]] = true;
665 roots[groups[j]] = i; //set root to self to start
667 //while you aren't at root
668 while(t->tree[index].getParent() != -1){
670 if (m->control_pressed) { return roots; }
672 //do both your chidren have have descendants from the users groups?
673 int lc = t->tree[index].getLChild();
674 int rc = t->tree[index].getRChild();
677 map<string, int>:: iterator itGroup = t->tree[lc].pcount.find(groups[j]);
678 if (itGroup != t->tree[lc].pcount.end()) { LpcountSize++; }
681 itGroup = t->tree[rc].pcount.find(groups[j]);
682 if (itGroup != t->tree[rc].pcount.end()) { RpcountSize++; }
684 if ((LpcountSize != 0) && (RpcountSize != 0)) { //possible root
685 roots[groups[j]] = index;
688 index = t->tree[index].getParent();
697 catch(exception& e) {
698 m->errorOut(e, "PhyloDiversityCommand", "getRootForGroups");
702 //**********************************************************************************************************************