2 * phylodiversitycommand.cpp
5 * Created by westcott on 4/30/10.
6 * Copyright 2010 Schloss Lab. All rights reserved.
10 #include "phylodiversitycommand.h"
12 //**********************************************************************************************************************
13 vector<string> PhyloDiversityCommand::setParameters(){
16 CommandParameter ptree("tree", "InputTypes", "", "", "none", "none", "none",false,true); parameters.push_back(ptree);
17 CommandParameter pgroup("group", "InputTypes", "", "", "none", "none", "none",false,true); parameters.push_back(pgroup);
18 CommandParameter pname("name", "InputTypes", "", "", "none", "none", "none",false,false); parameters.push_back(pname);
19 CommandParameter pgroups("groups", "String", "", "", "", "", "",false,false); parameters.push_back(pgroups);
20 CommandParameter piters("iters", "Number", "", "1000", "", "", "",false,false); parameters.push_back(piters);
21 CommandParameter pfreq("freq", "Number", "", "100", "", "", "",false,false); parameters.push_back(pfreq);
22 CommandParameter pprocessors("processors", "Number", "", "1", "", "", "",false,false); parameters.push_back(pprocessors);
23 CommandParameter prarefy("rarefy", "Boolean", "", "F", "", "", "",false,false); parameters.push_back(prarefy);
24 CommandParameter psummary("summary", "Boolean", "", "T", "", "", "",false,false); parameters.push_back(psummary);
25 CommandParameter pcollect("collect", "Boolean", "", "F", "", "", "",false,false); parameters.push_back(pcollect);
26 CommandParameter pscale("scale", "Boolean", "", "F", "", "", "",false,false); parameters.push_back(pscale);
27 CommandParameter pinputdir("inputdir", "String", "", "", "", "", "",false,false); parameters.push_back(pinputdir);
28 CommandParameter poutputdir("outputdir", "String", "", "", "", "", "",false,false); parameters.push_back(poutputdir);
30 vector<string> myArray;
31 for (int i = 0; i < parameters.size(); i++) { myArray.push_back(parameters[i].name); }
35 m->errorOut(e, "PhyloDiversityCommand", "setParameters");
39 //**********************************************************************************************************************
40 string PhyloDiversityCommand::getHelpString(){
42 string helpString = "";
43 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";
44 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";
45 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";
46 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";
47 helpString += "The scale parameter is used indicate that you want your output scaled to the number of sequences sampled, default = false. \n";
48 helpString += "The rarefy parameter allows you to create a rarefaction curve. The default is false.\n";
49 helpString += "The collect parameter allows you to create a collectors curve. The default is false.\n";
50 helpString += "The summary parameter allows you to create a .summary file. The default is true.\n";
51 helpString += "The processors parameter allows you to specify the number of processors to use. The default is 1.\n";
52 helpString += "The phylo.diversity command should be in the following format: phylo.diversity(groups=yourGroups, rarefy=yourRarefy, iters=yourIters).\n";
53 helpString += "Example phylo.diversity(groups=A-B-C, rarefy=T, iters=500).\n";
54 helpString += "The phylo.diversity command output two files: .phylo.diversity and if rarefy=T, .rarefaction.\n";
55 helpString += "Note: No spaces between parameter labels (i.e. groups), '=' and parameters (i.e.yourGroups).\n";
59 m->errorOut(e, "PhyloDiversityCommand", "getHelpString");
65 //**********************************************************************************************************************
66 PhyloDiversityCommand::PhyloDiversityCommand(){
68 abort = true; calledHelp = true;
70 vector<string> tempOutNames;
71 outputTypes["phylodiv"] = tempOutNames;
72 outputTypes["rarefy"] = tempOutNames;
73 outputTypes["summary"] = tempOutNames;
76 m->errorOut(e, "PhyloDiversityCommand", "PhyloDiversityCommand");
80 //**********************************************************************************************************************
81 PhyloDiversityCommand::PhyloDiversityCommand(string option) {
83 abort = false; calledHelp = false;
85 //allow user to run help
86 if(option == "help") { help(); abort = true; calledHelp = true; }
87 else if(option == "citation") { citation(); abort = true; calledHelp = true;}
90 vector<string> myArray = setParameters();;
92 OptionParser parser(option);
93 map<string,string> parameters = parser.getParameters();
94 map<string,string>::iterator it;
96 ValidParameters validParameter;
98 //check to make sure all parameters are valid for command
99 for (map<string,string>::iterator it = parameters.begin(); it != parameters.end(); it++) {
100 if (validParameter.isValidParameter(it->first, myArray, it->second) != true) { abort = true; }
103 //initialize outputTypes
104 vector<string> tempOutNames;
105 outputTypes["phylodiv"] = tempOutNames;
106 outputTypes["rarefy"] = tempOutNames;
107 outputTypes["summary"] = tempOutNames;
109 //if the user changes the input directory command factory will send this info to us in the output parameter
110 string inputDir = validParameter.validFile(parameters, "inputdir", false);
111 if (inputDir == "not found"){ inputDir = ""; }
114 it = parameters.find("tree");
115 //user has given a template file
116 if(it != parameters.end()){
117 path = m->hasPath(it->second);
118 //if the user has not given a path then, add inputdir. else leave path alone.
119 if (path == "") { parameters["tree"] = inputDir + it->second; }
122 it = parameters.find("group");
123 //user has given a template file
124 if(it != parameters.end()){
125 path = m->hasPath(it->second);
126 //if the user has not given a path then, add inputdir. else leave path alone.
127 if (path == "") { parameters["group"] = inputDir + it->second; }
130 it = parameters.find("name");
131 //user has given a template file
132 if(it != parameters.end()){
133 path = m->hasPath(it->second);
134 //if the user has not given a path then, add inputdir. else leave path alone.
135 if (path == "") { parameters["name"] = inputDir + it->second; }
141 m->namesOfGroups.clear();
142 m->Treenames.clear();
145 //check for required parameters
146 treefile = validParameter.validFile(parameters, "tree", true);
147 if (treefile == "not open") { abort = true; }
148 else if (treefile == "not found") {
149 //if there is a current design file, use it
150 treefile = m->getTreeFile();
151 if (treefile != "") { m->mothurOut("Using " + treefile + " as input file for the tree parameter."); m->mothurOutEndLine(); }
152 else { m->mothurOut("You have no current tree file and the tree parameter is required."); m->mothurOutEndLine(); abort = true; }
153 }else { m->setTreeFile(treefile); }
155 //check for required parameters
156 groupfile = validParameter.validFile(parameters, "group", true);
157 if (groupfile == "not open") { groupfile = ""; abort = true; }
158 else if (groupfile == "not found") { groupfile = ""; }
159 else { m->setGroupFile(groupfile); }
161 namefile = validParameter.validFile(parameters, "name", true);
162 if (namefile == "not open") { abort = true; }
163 else if (namefile == "not found") { namefile = ""; }
164 else { m->setNameFile(namefile); }
166 outputDir = validParameter.validFile(parameters, "outputdir", false); if (outputDir == "not found"){ outputDir = m->hasPath(treefile); }
169 temp = validParameter.validFile(parameters, "freq", false); if (temp == "not found") { temp = "100"; }
172 temp = validParameter.validFile(parameters, "iters", false); if (temp == "not found") { temp = "1000"; }
173 convert(temp, iters);
175 temp = validParameter.validFile(parameters, "rarefy", false); if (temp == "not found") { temp = "F"; }
176 rarefy = m->isTrue(temp);
177 if (!rarefy) { iters = 1; }
179 temp = validParameter.validFile(parameters, "summary", false); if (temp == "not found") { temp = "T"; }
180 summary = m->isTrue(temp);
182 temp = validParameter.validFile(parameters, "scale", false); if (temp == "not found") { temp = "F"; }
183 scale = m->isTrue(temp);
185 temp = validParameter.validFile(parameters, "collect", false); if (temp == "not found") { temp = "F"; }
186 collect = m->isTrue(temp);
188 temp = validParameter.validFile(parameters, "processors", false); if (temp == "not found"){ temp = m->getProcessors(); }
189 m->setProcessors(temp);
190 convert(temp, processors);
192 groups = validParameter.validFile(parameters, "groups", false);
193 if (groups == "not found") { groups = ""; }
195 m->splitAtDash(groups, Groups);
199 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; }
203 catch(exception& e) {
204 m->errorOut(e, "PhyloDiversityCommand", "PhyloDiversityCommand");
208 //**********************************************************************************************************************
210 int PhyloDiversityCommand::execute(){
213 if (abort == true) { if (calledHelp) { return 0; } return 2; }
215 m->setTreeFile(treefile);
217 if (groupfile != "") {
218 //read in group map info.
219 tmap = new TreeMap(groupfile);
221 }else{ //fake out by putting everyone in one group
222 Tree* tree = new Tree(treefile); delete tree; //extracts names from tree to make faked out groupmap
223 tmap = new TreeMap();
225 for (int i = 0; i < m->Treenames.size(); i++) { tmap->addSeq(m->Treenames[i], "Group1"); }
228 if (namefile != "") { readNamesFile(); }
230 read = new ReadNewickTree(treefile);
231 int readOk = read->read(tmap);
233 if (readOk != 0) { m->mothurOut("Read Terminated."); m->mothurOutEndLine(); delete tmap; delete read; return 0; }
235 read->AssembleTrees();
236 vector<Tree*> trees = read->getTrees();
239 //make sure all files match
240 //if you provide a namefile we will use the numNames in the namefile as long as the number of unique match the tree names size.
242 if (namefile != "") {
243 if (numUniquesInName == m->Treenames.size()) { numNamesInTree = nameMap.size(); }
244 else { numNamesInTree = m->Treenames.size(); }
245 }else { numNamesInTree = m->Treenames.size(); }
248 //output any names that are in group file but not in tree
249 if (numNamesInTree < tmap->getNumSeqs()) {
250 for (int i = 0; i < tmap->namesOfSeqs.size(); i++) {
251 //is that name in the tree?
253 for (int j = 0; j < m->Treenames.size(); j++) {
254 if (tmap->namesOfSeqs[i] == m->Treenames[j]) { break; } //found it
258 if (m->control_pressed) {
259 delete tmap; for (int i = 0; i < trees.size(); i++) { delete trees[i]; }
260 for (int i = 0; i < outputNames.size(); i++) { m->mothurRemove(outputNames[i]); } outputTypes.clear();
265 //then you did not find it so report it
266 if (count == m->Treenames.size()) {
267 //if it is in your namefile then don't remove
268 map<string, string>::iterator it = nameMap.find(tmap->namesOfSeqs[i]);
270 if (it == nameMap.end()) {
271 m->mothurOut(tmap->namesOfSeqs[i] + " is in your groupfile and not in your tree. It will be disregarded."); m->mothurOutEndLine();
272 tmap->removeSeq(tmap->namesOfSeqs[i]);
273 i--; //need this because removeSeq removes name from namesOfSeqs
279 SharedUtil* util = new SharedUtil();
280 util->setGroups(m->Groups, tmap->namesOfGroups, "phylo.diversity"); //sets the groups the user wants to analyze
283 //incase the user had some mismatches between the tree and group files we don't want group xxx to be analyzed
284 for (int i = 0; i < m->Groups.size(); i++) { if (m->Groups[i] == "xxx") { m->Groups.erase(m->Groups.begin()+i); break; } }
286 vector<string> outputNames;
288 //for each of the users trees
289 for(int i = 0; i < trees.size(); i++) {
291 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; }
293 ofstream outSum, outRare, outCollect;
294 string outSumFile = outputDir + m->getRootName(m->getSimpleName(treefile)) + toString(i+1) + ".phylodiv.summary";
295 string outRareFile = outputDir + m->getRootName(m->getSimpleName(treefile)) + toString(i+1) + ".phylodiv.rarefaction";
296 string outCollectFile = outputDir + m->getRootName(m->getSimpleName(treefile)) + toString(i+1) + ".phylodiv";
298 if (summary) { m->openOutputFile(outSumFile, outSum); outputNames.push_back(outSumFile); outputTypes["summary"].push_back(outSumFile); }
299 if (rarefy) { m->openOutputFile(outRareFile, outRare); outputNames.push_back(outRareFile); outputTypes["rarefy"].push_back(outRareFile); }
300 if (collect) { m->openOutputFile(outCollectFile, outCollect); outputNames.push_back(outCollectFile); outputTypes["phylodiv"].push_back(outCollectFile); }
302 int numLeafNodes = trees[i]->getNumLeaves();
304 //create a vector containing indexes of leaf nodes, randomize it, select nodes to send to calculator
305 vector<int> randomLeaf;
306 for (int j = 0; j < numLeafNodes; j++) {
307 if (m->inUsersGroups(trees[i]->tree[j].getGroup(), m->Groups) == true) { //is this a node from the group the user selected.
308 randomLeaf.push_back(j);
312 numLeafNodes = randomLeaf.size(); //reset the number of leaf nodes you are using
314 //each group, each sampling, if no rarefy iters = 1;
315 map<string, vector<float> > diversity;
317 //each group, each sampling, if no rarefy iters = 1;
318 map<string, vector<float> > sumDiversity;
320 //find largest group total
321 int largestGroup = 0;
322 for (int j = 0; j < m->Groups.size(); j++) {
323 if (tmap->seqsPerGroup[m->Groups[j]] > largestGroup) { largestGroup = tmap->seqsPerGroup[m->Groups[j]]; }
325 //initialize diversity
326 diversity[m->Groups[j]].resize(tmap->seqsPerGroup[m->Groups[j]]+1, 0.0); //numSampled
329 //initialize sumDiversity
330 sumDiversity[m->Groups[j]].resize(tmap->seqsPerGroup[m->Groups[j]]+1, 0.0);
333 //convert freq percentage to number
335 if (freq < 1.0) { increment = largestGroup * freq;
336 }else { increment = freq; }
338 //initialize sampling spots
339 set<int> numSampledList;
340 for(int k = 1; k <= largestGroup; k++){ if((k == 1) || (k % increment == 0)){ numSampledList.insert(k); } }
341 if(largestGroup % increment != 0){ numSampledList.insert(largestGroup); }
343 //add other groups ending points
344 for (int j = 0; j < m->Groups.size(); j++) {
345 if (numSampledList.count(diversity[m->Groups[j]].size()-1) == 0) { numSampledList.insert(diversity[m->Groups[j]].size()-1); }
348 #if defined (__APPLE__) || (__MACH__) || (linux) || (__linux)
350 driver(trees[i], diversity, sumDiversity, iters, increment, randomLeaf, numSampledList, outCollect, outSum, true);
353 vector<int> procIters;
355 int numItersPerProcessor = iters / processors;
357 //divide iters between processes
358 for (int h = 0; h < processors; h++) {
359 if(h == processors - 1){
360 numItersPerProcessor = iters - h * numItersPerProcessor;
362 procIters.push_back(numItersPerProcessor);
365 createProcesses(procIters, trees[i], diversity, sumDiversity, iters, increment, randomLeaf, numSampledList, outCollect, outSum);
367 }else{ //no need to paralellize if you dont want to rarefy
368 driver(trees[i], diversity, sumDiversity, iters, increment, randomLeaf, numSampledList, outCollect, outSum, true);
373 driver(trees[i], diversity, sumDiversity, iters, increment, randomLeaf, numSampledList, outCollect, outSum, true);
376 if (rarefy) { printData(numSampledList, sumDiversity, outRare, iters); }
380 if (m->control_pressed) { for (int i = 0; i < outputNames.size(); i++) { m->mothurRemove(outputNames[i]); } return 0; }
382 m->mothurOutEndLine();
383 m->mothurOut("Output File Names: "); m->mothurOutEndLine();
384 for (int i = 0; i < outputNames.size(); i++) { m->mothurOut(outputNames[i]); m->mothurOutEndLine(); }
385 m->mothurOutEndLine();
390 catch(exception& e) {
391 m->errorOut(e, "PhyloDiversityCommand", "execute");
395 //**********************************************************************************************************************
396 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){
398 #if defined (__APPLE__) || (__MACH__) || (linux) || (__linux)
401 vector<int> processIDS;
402 map< string, vector<float> >::iterator itSum;
404 //loop through and create all the processes you want
405 while (process != processors) {
409 processIDS.push_back(pid); //create map from line number to pid so you can append files in correct order later
412 driver(t, div, sumDiv, procIters[process], increment, randomLeaf, numSampledList, outCollect, outSum, false);
414 string outTemp = outputDir + toString(getpid()) + ".sumDiv.temp";
416 m->openOutputFile(outTemp, out);
418 //output the sumDIversity
419 for (itSum = sumDiv.begin(); itSum != sumDiv.end(); itSum++) {
420 out << itSum->first << '\t' << (itSum->second).size() << '\t';
421 for (int k = 0; k < (itSum->second).size(); k++) {
422 out << (itSum->second)[k] << '\t';
431 m->mothurOut("[ERROR]: unable to spawn the necessary processes."); m->mothurOutEndLine();
432 for (int i = 0; i < processIDS.size(); i++) { kill (processIDS[i], SIGINT); }
437 driver(t, div, sumDiv, procIters[0], increment, randomLeaf, numSampledList, outCollect, outSum, true);
439 //force parent to wait until all the processes are done
440 for (int i=0;i<(processors-1);i++) {
441 int temp = processIDS[i];
445 //get data created by processes
446 for (int i=0;i<(processors-1);i++) {
448 //input the sumDIversity
449 string inTemp = outputDir + toString(processIDS[i]) + ".sumDiv.temp";
451 m->openInputFile(inTemp, in);
453 //output the sumDIversity
454 for (int j = 0; j < sumDiv.size(); j++) {
458 in >> group >> size; m->gobble(in);
460 for (int k = 0; k < size; k++) {
464 sumDiv[group][k] += tempVal;
470 m->mothurRemove(inTemp);
478 catch(exception& e) {
479 m->errorOut(e, "PhyloDiversityCommand", "createProcesses");
483 //**********************************************************************************************************************
484 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){
486 int numLeafNodes = randomLeaf.size();
488 for (int l = 0; l < numIters; l++) {
489 random_shuffle(randomLeaf.begin(), randomLeaf.end());
492 map<string, int> counts;
493 map< string, set<int> > countedBranch;
494 for (int j = 0; j < m->Groups.size(); j++) { counts[m->Groups[j]] = 0; countedBranch[m->Groups[j]].insert(-2); } //add dummy index to initialize countedBranch sets
496 for(int k = 0; k < numLeafNodes; k++){
498 if (m->control_pressed) { return 0; }
500 //calc branch length of randomLeaf k
501 vector<float> br = calcBranchLength(t, randomLeaf[k], countedBranch);
503 //for each group in the groups update the total branch length accounting for the names file
504 vector<string> groups = t->tree[randomLeaf[k]].getGroup();
506 for (int j = 0; j < groups.size(); j++) {
507 int numSeqsInGroupJ = 0;
508 map<string, int>::iterator it;
509 it = t->tree[randomLeaf[k]].pcount.find(groups[j]);
510 if (it != t->tree[randomLeaf[k]].pcount.end()) { //this leaf node contains seqs from group j
511 numSeqsInGroupJ = it->second;
514 if (numSeqsInGroupJ != 0) { div[groups[j]][(counts[groups[j]]+1)] = div[groups[j]][counts[groups[j]]] + br[j]; }
516 for (int s = (counts[groups[j]]+2); s <= (counts[groups[j]]+numSeqsInGroupJ); s++) {
517 div[groups[j]][s] = div[groups[j]][s-1]; //update counts, but don't add in redundant branch lengths
519 counts[groups[j]] += numSeqsInGroupJ;
524 //add this diversity to the sum
525 for (int j = 0; j < m->Groups.size(); j++) {
526 for (int g = 0; g < div[m->Groups[j]].size(); g++) {
527 sumDiv[m->Groups[j]][g] += div[m->Groups[j]][g];
532 if ((collect) && (l == 0) && doSumCollect) { printData(numSampledList, div, outCollect, 1); }
533 if ((summary) && (l == 0) && doSumCollect) { printSumData(div, outSum, 1); }
539 catch(exception& e) {
540 m->errorOut(e, "PhyloDiversityCommand", "driver");
545 //**********************************************************************************************************************
547 void PhyloDiversityCommand::printSumData(map< string, vector<float> >& div, ofstream& out, int numIters){
550 out << "Groups\tnumSampled\tphyloDiversity" << endl;
552 out.setf(ios::fixed, ios::floatfield); out.setf(ios::showpoint);
554 for (int j = 0; j < m->Groups.size(); j++) {
555 int numSampled = (div[m->Groups[j]].size()-1);
556 out << m->Groups[j] << '\t' << numSampled << '\t';
560 if (scale) { score = (div[m->Groups[j]][numSampled] / (float)numIters) / (float)numSampled; }
561 else { score = div[m->Groups[j]][numSampled] / (float)numIters; }
563 out << setprecision(4) << score << endl;
569 catch(exception& e) {
570 m->errorOut(e, "PhyloDiversityCommand", "printSumData");
574 //**********************************************************************************************************************
576 void PhyloDiversityCommand::printData(set<int>& num, map< string, vector<float> >& div, ofstream& out, int numIters){
579 out << "numSampled\t";
580 for (int i = 0; i < m->Groups.size(); i++) { out << m->Groups[i] << '\t'; }
583 out.setf(ios::fixed, ios::floatfield); out.setf(ios::showpoint);
585 for (set<int>::iterator it = num.begin(); it != num.end(); it++) {
586 int numSampled = *it;
588 out << numSampled << '\t';
590 for (int j = 0; j < m->Groups.size(); j++) {
591 if (numSampled < div[m->Groups[j]].size()) {
593 if (scale) { score = (div[m->Groups[j]][numSampled] / (float)numIters) / (float)numSampled; }
594 else { score = div[m->Groups[j]][numSampled] / (float)numIters; }
596 out << setprecision(4) << score << '\t';
597 }else { out << "NA" << '\t'; }
605 catch(exception& e) {
606 m->errorOut(e, "PhyloDiversityCommand", "printData");
610 //**********************************************************************************************************************
611 //need a vector of floats one branch length for every group the node represents.
612 vector<float> PhyloDiversityCommand::calcBranchLength(Tree* t, int leaf, map< string, set<int> >& counted){
615 //calc the branch length
616 //while you aren't at root
620 vector<string> groups = t->tree[leaf].getGroup();
621 sums.resize(groups.size(), 0.0);
623 map<string, map<int, double> > tempTotals; //maps node to total Branch Length
624 map< string, set<int> > tempCounted;
625 set<int>::iterator it;
628 if(t->tree[index].getBranchLength() != -1){
629 for (int k = 0; k < groups.size(); k++) {
630 sums[k] += abs(t->tree[index].getBranchLength());
631 counted[groups[k]].insert(index);
635 for (int k = 0; k < groups.size(); k++) {
636 tempTotals[groups[k]][index] = 0.0;
639 index = t->tree[index].getParent();
641 //while you aren't at root
642 while(t->tree[index].getParent() != -1){
644 if (m->control_pressed) { return sums; }
647 for (int k = 0; k < groups.size(); k++) {
648 map<string, int>::iterator itGroup = t->tree[index].pcount.find(groups[k]);
649 if (itGroup != t->tree[index].pcount.end()) { pcountSize++; }
651 //do both your chidren have have descendants from the users groups?
652 int lc = t->tree[index].getLChild();
653 int rc = t->tree[index].getRChild();
656 itGroup = t->tree[lc].pcount.find(groups[k]);
657 if (itGroup != t->tree[lc].pcount.end()) { LpcountSize++; }
660 itGroup = t->tree[rc].pcount.find(groups[k]);
661 if (itGroup != t->tree[rc].pcount.end()) { RpcountSize++; }
663 //if yes, add your childrens tempTotals
664 if ((LpcountSize != 0) && (RpcountSize != 0)) {
665 sums[k] += tempTotals[groups[k]][lc] + tempTotals[groups[k]][rc];
667 for (it = tempCounted[groups[k]].begin(); it != tempCounted[groups[k]].end(); it++) { counted[groups[k]].insert(*it); }
669 //cout << "added to total " << tempTotals[lc] << '\t' << tempTotals[rc] << endl;
670 if (t->tree[index].getBranchLength() != -1) {
671 if (counted[groups[k]].count(index) == 0) {
672 tempTotals[groups[k]][index] = abs(t->tree[index].getBranchLength());
673 tempCounted[groups[k]].insert(index);
675 tempTotals[groups[k]][index] = 0.0;
678 tempTotals[groups[k]][index] = 0.0;
680 }else { //if no, your tempTotal is your childrens temp totals + your branch length
681 tempTotals[groups[k]][index] = tempTotals[groups[k]][lc] + tempTotals[groups[k]][rc];
683 if (counted[groups[k]].count(index) == 0) {
684 tempTotals[groups[k]][index] += abs(t->tree[index].getBranchLength());
685 tempCounted[groups[k]].insert(index);
689 //cout << "temptotal = "<< tempTotals[i] << endl;
692 index = t->tree[index].getParent();
698 catch(exception& e) {
699 m->errorOut(e, "PhyloDiversityCommand", "calcBranchLength");
703 /*****************************************************************/
704 int PhyloDiversityCommand::readNamesFile() {
707 numUniquesInName = 0;
710 m->openInputFile(namefile, in);
712 string first, second;
713 map<string, string>::iterator itNames;
716 in >> first >> second; m->gobble(in);
720 itNames = m->names.find(first);
721 if (itNames == m->names.end()) {
722 m->names[first] = second;
724 //we need a list of names in your namefile to use above when removing extra seqs above so we don't remove them
725 vector<string> dupNames;
726 m->splitAtComma(second, dupNames);
728 for (int i = 0; i < dupNames.size(); i++) {
729 nameMap[dupNames[i]] = dupNames[i];
730 if ((groupfile == "") && (i != 0)) { tmap->addSeq(dupNames[i], "Group1"); }
732 }else { m->mothurOut(first + " has already been seen in namefile, disregarding names file."); m->mothurOutEndLine(); in.close(); m->names.clear(); namefile = ""; return 1; }
738 catch(exception& e) {
739 m->errorOut(e, "PhyloDiversityCommand", "readNamesFile");
744 //**********************************************************************************************************************