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 pname("name", "InputTypes", "", "", "NameCount", "none", "none",false,false); parameters.push_back(pname);
19 CommandParameter pcount("count", "InputTypes", "", "", "NameCount-CountGroup", "none", "none",false,false); parameters.push_back(pcount);
20 CommandParameter pgroup("group", "InputTypes", "", "", "CountGroup", "none", "none",false,false); parameters.push_back(pgroup);
21 CommandParameter pgroups("groups", "String", "", "", "", "", "",false,false); parameters.push_back(pgroups);
22 CommandParameter piters("iters", "Number", "", "1000", "", "", "",false,false); parameters.push_back(piters);
23 CommandParameter pfreq("freq", "Number", "", "100", "", "", "",false,false); parameters.push_back(pfreq);
24 CommandParameter pprocessors("processors", "Number", "", "1", "", "", "",false,false); parameters.push_back(pprocessors);
25 CommandParameter prarefy("rarefy", "Boolean", "", "F", "", "", "",false,false); parameters.push_back(prarefy);
26 CommandParameter psummary("summary", "Boolean", "", "T", "", "", "",false,false); parameters.push_back(psummary);
27 CommandParameter pcollect("collect", "Boolean", "", "F", "", "", "",false,false); parameters.push_back(pcollect);
28 CommandParameter pscale("scale", "Boolean", "", "F", "", "", "",false,false); parameters.push_back(pscale);
29 CommandParameter pinputdir("inputdir", "String", "", "", "", "", "",false,false); parameters.push_back(pinputdir);
30 CommandParameter poutputdir("outputdir", "String", "", "", "", "", "",false,false); parameters.push_back(poutputdir);
32 vector<string> myArray;
33 for (int i = 0; i < parameters.size(); i++) { myArray.push_back(parameters[i].name); }
37 m->errorOut(e, "PhyloDiversityCommand", "setParameters");
41 //**********************************************************************************************************************
42 string PhyloDiversityCommand::getHelpString(){
44 string helpString = "";
45 helpString += "The phylo.diversity command parameters are tree, group, name, count, groups, iters, freq, processors, scale, rarefy, collect and summary. tree and group are required, unless you have valid current files.\n";
46 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";
47 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";
48 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";
49 helpString += "The scale parameter is used indicate that you want your output scaled to the number of sequences sampled, default = false. \n";
50 helpString += "The rarefy parameter allows you to create a rarefaction curve. The default is false.\n";
51 helpString += "The collect parameter allows you to create a collectors curve. The default is false.\n";
52 helpString += "The summary parameter allows you to create a .summary file. The default is true.\n";
53 helpString += "The processors parameter allows you to specify the number of processors to use. The default is 1.\n";
54 helpString += "The phylo.diversity command should be in the following format: phylo.diversity(groups=yourGroups, rarefy=yourRarefy, iters=yourIters).\n";
55 helpString += "Example phylo.diversity(groups=A-B-C, rarefy=T, iters=500).\n";
56 helpString += "The phylo.diversity command output two files: .phylo.diversity and if rarefy=T, .rarefaction.\n";
57 helpString += "Note: No spaces between parameter labels (i.e. groups), '=' and parameters (i.e.yourGroups).\n";
61 m->errorOut(e, "PhyloDiversityCommand", "getHelpString");
65 //**********************************************************************************************************************
66 string PhyloDiversityCommand::getOutputFileNameTag(string type, string inputName=""){
68 string outputFileName = "";
69 map<string, vector<string> >::iterator it;
71 //is this a type this command creates
72 it = outputTypes.find(type);
73 if (it == outputTypes.end()) { m->mothurOut("[ERROR]: this command doesn't create a " + type + " output file.\n"); }
75 if (type == "phylodiv") { outputFileName = "phylodiv"; }
76 else if (type == "rarefy") { outputFileName = "phylodiv.rarefaction"; }
77 else if (type == "summary") { outputFileName = "phylodiv.summary"; }
78 else { m->mothurOut("[ERROR]: No definition for type " + type + " output file tag.\n"); m->control_pressed = true; }
80 return outputFileName;
83 m->errorOut(e, "PhyloDiversityCommand", "getOutputFileNameTag");
88 //**********************************************************************************************************************
89 PhyloDiversityCommand::PhyloDiversityCommand(){
91 abort = true; calledHelp = true;
93 vector<string> tempOutNames;
94 outputTypes["phylodiv"] = tempOutNames;
95 outputTypes["rarefy"] = tempOutNames;
96 outputTypes["summary"] = tempOutNames;
99 m->errorOut(e, "PhyloDiversityCommand", "PhyloDiversityCommand");
103 //**********************************************************************************************************************
104 PhyloDiversityCommand::PhyloDiversityCommand(string option) {
106 abort = false; calledHelp = false;
108 //allow user to run help
109 if(option == "help") { help(); abort = true; calledHelp = true; }
110 else if(option == "citation") { citation(); abort = true; calledHelp = true;}
113 vector<string> myArray = setParameters();;
115 OptionParser parser(option);
116 map<string,string> parameters = parser.getParameters();
117 map<string,string>::iterator it;
119 ValidParameters validParameter;
121 //check to make sure all parameters are valid for command
122 for (map<string,string>::iterator it = parameters.begin(); it != parameters.end(); it++) {
123 if (validParameter.isValidParameter(it->first, myArray, it->second) != true) { abort = true; }
126 //initialize outputTypes
127 vector<string> tempOutNames;
128 outputTypes["phylodiv"] = tempOutNames;
129 outputTypes["rarefy"] = tempOutNames;
130 outputTypes["summary"] = tempOutNames;
132 //if the user changes the input directory command factory will send this info to us in the output parameter
133 string inputDir = validParameter.validFile(parameters, "inputdir", false);
134 if (inputDir == "not found"){ inputDir = ""; }
137 it = parameters.find("tree");
138 //user has given a template file
139 if(it != parameters.end()){
140 path = m->hasPath(it->second);
141 //if the user has not given a path then, add inputdir. else leave path alone.
142 if (path == "") { parameters["tree"] = inputDir + it->second; }
145 it = parameters.find("group");
146 //user has given a template file
147 if(it != parameters.end()){
148 path = m->hasPath(it->second);
149 //if the user has not given a path then, add inputdir. else leave path alone.
150 if (path == "") { parameters["group"] = inputDir + it->second; }
153 it = parameters.find("name");
154 //user has given a template file
155 if(it != parameters.end()){
156 path = m->hasPath(it->second);
157 //if the user has not given a path then, add inputdir. else leave path alone.
158 if (path == "") { parameters["name"] = inputDir + it->second; }
161 it = parameters.find("count");
162 //user has given a template file
163 if(it != parameters.end()){
164 path = m->hasPath(it->second);
165 //if the user has not given a path then, add inputdir. else leave path alone.
166 if (path == "") { parameters["count"] = inputDir + it->second; }
170 //check for required parameters
171 treefile = validParameter.validFile(parameters, "tree", true);
172 if (treefile == "not open") { treefile = ""; abort = true; }
173 else if (treefile == "not found") {
174 //if there is a current design file, use it
175 treefile = m->getTreeFile();
176 if (treefile != "") { m->mothurOut("Using " + treefile + " as input file for the tree parameter."); m->mothurOutEndLine(); }
177 else { m->mothurOut("You have no current tree file and the tree parameter is required."); m->mothurOutEndLine(); abort = true; }
178 }else { m->setTreeFile(treefile); }
180 //check for required parameters
181 groupfile = validParameter.validFile(parameters, "group", true);
182 if (groupfile == "not open") { groupfile = ""; abort = true; }
183 else if (groupfile == "not found") { groupfile = ""; }
184 else { m->setGroupFile(groupfile); }
186 namefile = validParameter.validFile(parameters, "name", true);
187 if (namefile == "not open") { namefile = ""; abort = true; }
188 else if (namefile == "not found") { namefile = ""; }
189 else { m->setNameFile(namefile); }
191 countfile = validParameter.validFile(parameters, "count", true);
192 if (countfile == "not open") { countfile = ""; abort = true; }
193 else if (countfile == "not found") { countfile = ""; }
194 else { m->setCountTableFile(countfile); }
196 if ((namefile != "") && (countfile != "")) {
197 m->mothurOut("[ERROR]: you may only use one of the following: name or count."); m->mothurOutEndLine(); abort = true;
200 if ((groupfile != "") && (countfile != "")) {
201 m->mothurOut("[ERROR]: you may only use one of the following: group or count."); m->mothurOutEndLine(); abort=true;
204 outputDir = validParameter.validFile(parameters, "outputdir", false); if (outputDir == "not found"){ outputDir = m->hasPath(treefile); }
207 temp = validParameter.validFile(parameters, "freq", false); if (temp == "not found") { temp = "100"; }
208 m->mothurConvert(temp, freq);
210 temp = validParameter.validFile(parameters, "iters", false); if (temp == "not found") { temp = "1000"; }
211 m->mothurConvert(temp, iters);
213 temp = validParameter.validFile(parameters, "rarefy", false); if (temp == "not found") { temp = "F"; }
214 rarefy = m->isTrue(temp);
215 if (!rarefy) { iters = 1; }
217 temp = validParameter.validFile(parameters, "summary", false); if (temp == "not found") { temp = "T"; }
218 summary = m->isTrue(temp);
220 temp = validParameter.validFile(parameters, "scale", false); if (temp == "not found") { temp = "F"; }
221 scale = m->isTrue(temp);
223 temp = validParameter.validFile(parameters, "collect", false); if (temp == "not found") { temp = "F"; }
224 collect = m->isTrue(temp);
226 temp = validParameter.validFile(parameters, "processors", false); if (temp == "not found"){ temp = m->getProcessors(); }
227 m->setProcessors(temp);
228 m->mothurConvert(temp, processors);
230 groups = validParameter.validFile(parameters, "groups", false);
231 if (groups == "not found") { groups = ""; }
233 m->splitAtDash(groups, Groups);
234 m->setGroups(Groups);
237 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; }
240 if (namefile == "") {
241 vector<string> files; files.push_back(treefile);
242 parser.getNameFile(files);
248 catch(exception& e) {
249 m->errorOut(e, "PhyloDiversityCommand", "PhyloDiversityCommand");
253 //**********************************************************************************************************************
255 int PhyloDiversityCommand::execute(){
258 if (abort == true) { if (calledHelp) { return 0; } return 2; }
260 int start = time(NULL);
262 m->setTreeFile(treefile);
264 if (countfile == "") { reader = new TreeReader(treefile, groupfile, namefile); }
265 else { reader = new TreeReader(treefile, countfile); }
266 vector<Tree*> trees = reader->getTrees();
267 ct = trees[0]->getCountTable();
271 vector<string> mGroups = m->getGroups();
272 vector<string> tGroups = ct->getNamesOfGroups();
273 util.setGroups(mGroups, tGroups, "phylo.diversity"); //sets the groups the user wants to analyze
275 //incase the user had some mismatches between the tree and group files we don't want group xxx to be analyzed
276 for (int i = 0; i < mGroups.size(); i++) { if (mGroups[i] == "xxx") { mGroups.erase(mGroups.begin()+i); break; } }
277 m->setGroups(mGroups);
279 vector<string> outputNames;
281 //for each of the users trees
282 for(int i = 0; i < trees.size(); i++) {
284 if (m->control_pressed) { delete ct; 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; }
286 ofstream outSum, outRare, outCollect;
287 string outSumFile = outputDir + m->getRootName(m->getSimpleName(treefile)) + toString(i+1) + "." + getOutputFileNameTag("summary");
288 string outRareFile = outputDir + m->getRootName(m->getSimpleName(treefile)) + toString(i+1) + "." + getOutputFileNameTag("rarefy");
289 string outCollectFile = outputDir + m->getRootName(m->getSimpleName(treefile)) + toString(i+1) + "." + getOutputFileNameTag("phylodiv");
291 if (summary) { m->openOutputFile(outSumFile, outSum); outputNames.push_back(outSumFile); outputTypes["summary"].push_back(outSumFile); }
292 if (rarefy) { m->openOutputFile(outRareFile, outRare); outputNames.push_back(outRareFile); outputTypes["rarefy"].push_back(outRareFile); }
293 if (collect) { m->openOutputFile(outCollectFile, outCollect); outputNames.push_back(outCollectFile); outputTypes["phylodiv"].push_back(outCollectFile); }
295 int numLeafNodes = trees[i]->getNumLeaves();
297 //create a vector containing indexes of leaf nodes, randomize it, select nodes to send to calculator
298 vector<int> randomLeaf;
299 for (int j = 0; j < numLeafNodes; j++) {
300 if (m->inUsersGroups(trees[i]->tree[j].getGroup(), mGroups) == true) { //is this a node from the group the user selected.
301 randomLeaf.push_back(j);
305 numLeafNodes = randomLeaf.size(); //reset the number of leaf nodes you are using
307 //each group, each sampling, if no rarefy iters = 1;
308 map<string, vector<float> > diversity;
310 //each group, each sampling, if no rarefy iters = 1;
311 map<string, vector<float> > sumDiversity;
313 //find largest group total
314 int largestGroup = 0;
315 for (int j = 0; j < mGroups.size(); j++) {
316 int numSeqsThisGroup = ct->getGroupCount(mGroups[j]);
317 if (numSeqsThisGroup > largestGroup) { largestGroup = numSeqsThisGroup; }
319 //initialize diversity
320 diversity[mGroups[j]].resize(numSeqsThisGroup+1, 0.0); //numSampled
323 //initialize sumDiversity
324 sumDiversity[mGroups[j]].resize(numSeqsThisGroup+1, 0.0);
327 //convert freq percentage to number
329 if (freq < 1.0) { increment = largestGroup * freq;
330 }else { increment = freq; }
332 //initialize sampling spots
333 set<int> numSampledList;
334 for(int k = 1; k <= largestGroup; k++){ if((k == 1) || (k % increment == 0)){ numSampledList.insert(k); } }
335 if(largestGroup % increment != 0){ numSampledList.insert(largestGroup); }
337 //add other groups ending points
338 for (int j = 0; j < mGroups.size(); j++) {
339 if (numSampledList.count(diversity[mGroups[j]].size()-1) == 0) { numSampledList.insert(diversity[mGroups[j]].size()-1); }
342 #if defined (__APPLE__) || (__MACH__) || (linux) || (__linux) || (__linux__) || (__unix__) || (__unix)
344 driver(trees[i], diversity, sumDiversity, iters, increment, randomLeaf, numSampledList, outCollect, outSum, true);
347 vector<int> procIters;
349 int numItersPerProcessor = iters / processors;
351 //divide iters between processes
352 for (int h = 0; h < processors; h++) {
353 if(h == processors - 1){
354 numItersPerProcessor = iters - h * numItersPerProcessor;
356 procIters.push_back(numItersPerProcessor);
359 createProcesses(procIters, trees[i], diversity, sumDiversity, iters, increment, randomLeaf, numSampledList, outCollect, outSum);
361 }else{ //no need to paralellize if you dont want to rarefy
362 driver(trees[i], diversity, sumDiversity, iters, increment, randomLeaf, numSampledList, outCollect, outSum, true);
367 driver(trees[i], diversity, sumDiversity, iters, increment, randomLeaf, numSampledList, outCollect, outSum, true);
370 if (rarefy) { printData(numSampledList, sumDiversity, outRare, iters); }
374 if (m->control_pressed) { for (int i = 0; i < outputNames.size(); i++) { m->mothurRemove(outputNames[i]); } return 0; }
376 m->mothurOut("It took " + toString(time(NULL) - start) + " secs to run phylo.diversity."); m->mothurOutEndLine();
379 m->mothurOutEndLine();
380 m->mothurOut("Output File Names: "); m->mothurOutEndLine();
381 for (int i = 0; i < outputNames.size(); i++) { m->mothurOut(outputNames[i]); m->mothurOutEndLine(); }
382 m->mothurOutEndLine();
387 catch(exception& e) {
388 m->errorOut(e, "PhyloDiversityCommand", "execute");
392 //**********************************************************************************************************************
393 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){
395 #if defined (__APPLE__) || (__MACH__) || (linux) || (__linux) || (__linux__) || (__unix__) || (__unix)
398 vector<int> processIDS;
399 map< string, vector<float> >::iterator itSum;
401 //loop through and create all the processes you want
402 while (process != processors) {
406 processIDS.push_back(pid); //create map from line number to pid so you can append files in correct order later
409 driver(t, div, sumDiv, procIters[process], increment, randomLeaf, numSampledList, outCollect, outSum, false);
411 string outTemp = outputDir + toString(getpid()) + ".sumDiv.temp";
413 m->openOutputFile(outTemp, out);
415 //output the sumDIversity
416 for (itSum = sumDiv.begin(); itSum != sumDiv.end(); itSum++) {
417 out << itSum->first << '\t' << (itSum->second).size() << '\t';
418 for (int k = 0; k < (itSum->second).size(); k++) {
419 out << (itSum->second)[k] << '\t';
428 m->mothurOut("[ERROR]: unable to spawn the necessary processes."); m->mothurOutEndLine();
429 for (int i = 0; i < processIDS.size(); i++) { kill (processIDS[i], SIGINT); }
434 driver(t, div, sumDiv, procIters[0], increment, randomLeaf, numSampledList, outCollect, outSum, true);
436 //force parent to wait until all the processes are done
437 for (int i=0;i<(processors-1);i++) {
438 int temp = processIDS[i];
442 //get data created by processes
443 for (int i=0;i<(processors-1);i++) {
445 //input the sumDIversity
446 string inTemp = outputDir + toString(processIDS[i]) + ".sumDiv.temp";
448 m->openInputFile(inTemp, in);
450 //output the sumDIversity
451 for (int j = 0; j < sumDiv.size(); j++) {
455 in >> group >> size; m->gobble(in);
457 for (int k = 0; k < size; k++) {
461 sumDiv[group][k] += tempVal;
467 m->mothurRemove(inTemp);
475 catch(exception& e) {
476 m->errorOut(e, "PhyloDiversityCommand", "createProcesses");
480 //**********************************************************************************************************************
481 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){
483 int numLeafNodes = randomLeaf.size();
484 vector<string> mGroups = m->getGroups();
486 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.
488 for (int l = 0; l < numIters; l++) {
489 random_shuffle(randomLeaf.begin(), randomLeaf.end());
492 map<string, int> counts;
493 vector< map<string, bool> > countedBranch;
494 for (int i = 0; i < t->getNumNodes(); i++) {
495 map<string, bool> temp;
496 for (int j = 0; j < mGroups.size(); j++) { temp[mGroups[j]] = false; }
497 countedBranch.push_back(temp);
500 for (int j = 0; j < mGroups.size(); j++) { counts[mGroups[j]] = 0; }
502 for(int k = 0; k < numLeafNodes; k++){
504 if (m->control_pressed) { return 0; }
506 //calc branch length of randomLeaf k
507 vector<float> br = calcBranchLength(t, randomLeaf[k], countedBranch, rootForGroup);
509 //for each group in the groups update the total branch length accounting for the names file
510 vector<string> groups = t->tree[randomLeaf[k]].getGroup();
512 for (int j = 0; j < groups.size(); j++) {
514 if (m->inUsersGroups(groups[j], mGroups)) {
515 int numSeqsInGroupJ = 0;
516 map<string, int>::iterator it;
517 it = t->tree[randomLeaf[k]].pcount.find(groups[j]);
518 if (it != t->tree[randomLeaf[k]].pcount.end()) { //this leaf node contains seqs from group j
519 numSeqsInGroupJ = it->second;
522 if (numSeqsInGroupJ != 0) { div[groups[j]][(counts[groups[j]]+1)] = div[groups[j]][counts[groups[j]]] + br[j]; }
524 for (int s = (counts[groups[j]]+2); s <= (counts[groups[j]]+numSeqsInGroupJ); s++) {
525 div[groups[j]][s] = div[groups[j]][s-1]; //update counts, but don't add in redundant branch lengths
527 counts[groups[j]] += numSeqsInGroupJ;
533 //add this diversity to the sum
534 for (int j = 0; j < mGroups.size(); j++) {
535 for (int g = 0; g < div[mGroups[j]].size(); g++) {
536 sumDiv[mGroups[j]][g] += div[mGroups[j]][g];
541 if ((collect) && (l == 0) && doSumCollect) { printData(numSampledList, div, outCollect, 1); }
542 if ((summary) && (l == 0) && doSumCollect) { printSumData(div, outSum, 1); }
548 catch(exception& e) {
549 m->errorOut(e, "PhyloDiversityCommand", "driver");
554 //**********************************************************************************************************************
556 void PhyloDiversityCommand::printSumData(map< string, vector<float> >& div, ofstream& out, int numIters){
559 out << "Groups\tnumSampled\tphyloDiversity" << endl;
561 out.setf(ios::fixed, ios::floatfield); out.setf(ios::showpoint);
563 vector<string> mGroups = m->getGroups();
564 for (int j = 0; j < mGroups.size(); j++) {
565 int numSampled = (div[mGroups[j]].size()-1);
566 out << mGroups[j] << '\t' << numSampled << '\t';
570 if (scale) { score = (div[mGroups[j]][numSampled] / (float)numIters) / (float)numSampled; }
571 else { score = div[mGroups[j]][numSampled] / (float)numIters; }
573 out << setprecision(4) << score << endl;
579 catch(exception& e) {
580 m->errorOut(e, "PhyloDiversityCommand", "printSumData");
584 //**********************************************************************************************************************
586 void PhyloDiversityCommand::printData(set<int>& num, map< string, vector<float> >& div, ofstream& out, int numIters){
589 out << "numSampled\t";
590 vector<string> mGroups = m->getGroups();
591 for (int i = 0; i < mGroups.size(); i++) { out << mGroups[i] << '\t'; }
594 out.setf(ios::fixed, ios::floatfield); out.setf(ios::showpoint);
596 for (set<int>::iterator it = num.begin(); it != num.end(); it++) {
597 int numSampled = *it;
599 out << numSampled << '\t';
601 for (int j = 0; j < mGroups.size(); j++) {
602 if (numSampled < div[mGroups[j]].size()) {
604 if (scale) { score = (div[mGroups[j]][numSampled] / (float)numIters) / (float)numSampled; }
605 else { score = div[mGroups[j]][numSampled] / (float)numIters; }
607 out << setprecision(4) << score << '\t';
608 }else { out << "NA" << '\t'; }
616 catch(exception& e) {
617 m->errorOut(e, "PhyloDiversityCommand", "printData");
621 //**********************************************************************************************************************
622 //need a vector of floats one branch length for every group the node represents.
623 vector<float> PhyloDiversityCommand::calcBranchLength(Tree* t, int leaf, vector< map<string, bool> >& counted, map<string, int> roots){
626 //calc the branch length
627 //while you aren't at root
631 vector<string> groups = t->tree[leaf].getGroup();
632 sums.resize(groups.size(), 0.0);
636 if(t->tree[index].getBranchLength() != -1){
637 for (int k = 0; k < groups.size(); k++) {
638 sums[k] += abs(t->tree[index].getBranchLength());
643 index = t->tree[index].getParent();
645 //while you aren't at root
646 while(t->tree[index].getParent() != -1){
648 if (m->control_pressed) { return sums; }
650 for (int k = 0; k < groups.size(); k++) {
652 if (index >= roots[groups[k]]) { counted[index][groups[k]] = true; } //if you are at this groups "root", then say we are done
654 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
655 if (t->tree[index].getBranchLength() != -1) {
656 sums[k] += abs(t->tree[index].getBranchLength());
658 counted[index][groups[k]] = true;
661 index = t->tree[index].getParent();
667 catch(exception& e) {
668 m->errorOut(e, "PhyloDiversityCommand", "calcBranchLength");
672 //**********************************************************************************************************************
673 map<string, int> PhyloDiversityCommand::getRootForGroups(Tree* t){
675 map<string, int> roots; //maps group to root for group, may not be root of tree
676 map<string, bool> done;
678 //initialize root for all groups to -1
679 for (int k = 0; k < (t->getCountTable())->getNamesOfGroups().size(); k++) { done[(t->getCountTable())->getNamesOfGroups()[k]] = false; }
681 for (int i = 0; i < t->getNumLeaves(); i++) {
683 vector<string> groups = t->tree[i].getGroup();
685 int index = t->tree[i].getParent();
687 for (int j = 0; j < groups.size(); j++) {
689 if (done[groups[j]] == false) { //we haven't found the root for this group yet
691 done[groups[j]] = true;
692 roots[groups[j]] = i; //set root to self to start
694 //while you aren't at root
695 while(t->tree[index].getParent() != -1){
697 if (m->control_pressed) { return roots; }
699 //do both your chidren have have descendants from the users groups?
700 int lc = t->tree[index].getLChild();
701 int rc = t->tree[index].getRChild();
704 map<string, int>:: iterator itGroup = t->tree[lc].pcount.find(groups[j]);
705 if (itGroup != t->tree[lc].pcount.end()) { LpcountSize++; }
708 itGroup = t->tree[rc].pcount.find(groups[j]);
709 if (itGroup != t->tree[rc].pcount.end()) { RpcountSize++; }
711 if ((LpcountSize != 0) && (RpcountSize != 0)) { //possible root
712 roots[groups[j]] = index;
715 index = t->tree[index].getParent();
724 catch(exception& e) {
725 m->errorOut(e, "PhyloDiversityCommand", "getRootForGroups");
729 //**********************************************************************************************************************