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; }
89 vector<string> myArray = setParameters();;
91 OptionParser parser(option);
92 map<string,string> parameters = parser.getParameters();
93 map<string,string>::iterator it;
95 ValidParameters validParameter;
97 //check to make sure all parameters are valid for command
98 for (map<string,string>::iterator it = parameters.begin(); it != parameters.end(); it++) {
99 if (validParameter.isValidParameter(it->first, myArray, it->second) != true) { abort = true; }
102 //initialize outputTypes
103 vector<string> tempOutNames;
104 outputTypes["phylodiv"] = tempOutNames;
105 outputTypes["rarefy"] = tempOutNames;
106 outputTypes["summary"] = tempOutNames;
108 //if the user changes the input directory command factory will send this info to us in the output parameter
109 string inputDir = validParameter.validFile(parameters, "inputdir", false);
110 if (inputDir == "not found"){ inputDir = ""; }
113 it = parameters.find("tree");
114 //user has given a template file
115 if(it != parameters.end()){
116 path = m->hasPath(it->second);
117 //if the user has not given a path then, add inputdir. else leave path alone.
118 if (path == "") { parameters["tree"] = inputDir + it->second; }
121 it = parameters.find("group");
122 //user has given a template file
123 if(it != parameters.end()){
124 path = m->hasPath(it->second);
125 //if the user has not given a path then, add inputdir. else leave path alone.
126 if (path == "") { parameters["group"] = inputDir + it->second; }
129 it = parameters.find("name");
130 //user has given a template file
131 if(it != parameters.end()){
132 path = m->hasPath(it->second);
133 //if the user has not given a path then, add inputdir. else leave path alone.
134 if (path == "") { parameters["name"] = inputDir + it->second; }
140 m->namesOfGroups.clear();
141 m->Treenames.clear();
144 //check for required parameters
145 treefile = validParameter.validFile(parameters, "tree", true);
146 if (treefile == "not open") { abort = true; }
147 else if (treefile == "not found") {
148 //if there is a current design file, use it
149 treefile = m->getTreeFile();
150 if (treefile != "") { m->mothurOut("Using " + treefile + " as input file for the tree parameter."); m->mothurOutEndLine(); }
151 else { m->mothurOut("You have no current tree file and the tree parameter is required."); m->mothurOutEndLine(); abort = true; }
154 //check for required parameters
155 groupfile = validParameter.validFile(parameters, "group", true);
156 if (groupfile == "not open") { abort = true; }
157 else if (groupfile == "not found") {
158 //if there is a current design file, use it
159 groupfile = m->getGroupFile();
160 if (groupfile != "") { m->mothurOut("Using " + groupfile + " as input file for the group parameter."); m->mothurOutEndLine(); }
161 else { m->mothurOut("You have no current group file and the group parameter is required."); m->mothurOutEndLine(); abort = true; }
164 namefile = validParameter.validFile(parameters, "name", true);
165 if (namefile == "not open") { abort = true; }
166 else if (namefile == "not found") { namefile = ""; }
168 outputDir = validParameter.validFile(parameters, "outputdir", false); if (outputDir == "not found"){ outputDir = m->hasPath(treefile); }
171 temp = validParameter.validFile(parameters, "freq", false); if (temp == "not found") { temp = "100"; }
174 temp = validParameter.validFile(parameters, "iters", false); if (temp == "not found") { temp = "1000"; }
175 convert(temp, iters);
177 temp = validParameter.validFile(parameters, "rarefy", false); if (temp == "not found") { temp = "F"; }
178 rarefy = m->isTrue(temp);
179 if (!rarefy) { iters = 1; }
181 temp = validParameter.validFile(parameters, "summary", false); if (temp == "not found") { temp = "T"; }
182 summary = m->isTrue(temp);
184 temp = validParameter.validFile(parameters, "scale", false); if (temp == "not found") { temp = "F"; }
185 scale = m->isTrue(temp);
187 temp = validParameter.validFile(parameters, "collect", false); if (temp == "not found") { temp = "F"; }
188 collect = m->isTrue(temp);
190 temp = validParameter.validFile(parameters, "processors", false); if (temp == "not found"){ temp = m->getProcessors(); }
191 m->setProcessors(temp);
192 convert(temp, processors);
194 groups = validParameter.validFile(parameters, "groups", false);
195 if (groups == "not found") { groups = ""; }
197 m->splitAtDash(groups, Groups);
201 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; }
205 catch(exception& e) {
206 m->errorOut(e, "PhyloDiversityCommand", "PhyloDiversityCommand");
210 //**********************************************************************************************************************
212 int PhyloDiversityCommand::execute(){
215 if (abort == true) { if (calledHelp) { return 0; } return 2; }
217 m->setTreeFile(treefile);
219 //read in group map info.
220 tmap = new TreeMap(groupfile);
223 if (namefile != "") { readNamesFile(); }
225 read = new ReadNewickTree(treefile);
226 int readOk = read->read(tmap);
228 if (readOk != 0) { m->mothurOut("Read Terminated."); m->mothurOutEndLine(); delete tmap; delete read; return 0; }
230 read->AssembleTrees();
231 vector<Tree*> trees = read->getTrees();
234 //make sure all files match
235 //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.
237 if (namefile != "") {
238 if (numUniquesInName == m->Treenames.size()) { numNamesInTree = nameMap.size(); }
239 else { numNamesInTree = m->Treenames.size(); }
240 }else { numNamesInTree = m->Treenames.size(); }
243 //output any names that are in group file but not in tree
244 if (numNamesInTree < tmap->getNumSeqs()) {
245 for (int i = 0; i < tmap->namesOfSeqs.size(); i++) {
246 //is that name in the tree?
248 for (int j = 0; j < m->Treenames.size(); j++) {
249 if (tmap->namesOfSeqs[i] == m->Treenames[j]) { break; } //found it
253 if (m->control_pressed) {
254 delete tmap; for (int i = 0; i < trees.size(); i++) { delete trees[i]; }
255 for (int i = 0; i < outputNames.size(); i++) { remove(outputNames[i].c_str()); } outputTypes.clear();
260 //then you did not find it so report it
261 if (count == m->Treenames.size()) {
262 //if it is in your namefile then don't remove
263 map<string, string>::iterator it = nameMap.find(tmap->namesOfSeqs[i]);
265 if (it == nameMap.end()) {
266 m->mothurOut(tmap->namesOfSeqs[i] + " is in your groupfile and not in your tree. It will be disregarded."); m->mothurOutEndLine();
267 tmap->removeSeq(tmap->namesOfSeqs[i]);
268 i--; //need this because removeSeq removes name from namesOfSeqs
274 SharedUtil* util = new SharedUtil();
275 util->setGroups(m->Groups, tmap->namesOfGroups, "treegroup"); //sets the groups the user wants to analyze
278 //incase the user had some mismatches between the tree and group files we don't want group xxx to be analyzed
279 for (int i = 0; i < m->Groups.size(); i++) { if (m->Groups[i] == "xxx") { m->Groups.erase(m->Groups.begin()+i); break; } }
281 vector<string> outputNames;
283 //for each of the users trees
284 for(int i = 0; i < trees.size(); i++) {
286 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++) { remove(outputNames[j].c_str()); } return 0; }
288 ofstream outSum, outRare, outCollect;
289 string outSumFile = outputDir + m->getRootName(m->getSimpleName(treefile)) + toString(i+1) + ".phylodiv.summary";
290 string outRareFile = outputDir + m->getRootName(m->getSimpleName(treefile)) + toString(i+1) + ".phylodiv.rarefaction";
291 string outCollectFile = outputDir + m->getRootName(m->getSimpleName(treefile)) + toString(i+1) + ".phylodiv";
293 if (summary) { m->openOutputFile(outSumFile, outSum); outputNames.push_back(outSumFile); outputTypes["summary"].push_back(outSumFile); }
294 if (rarefy) { m->openOutputFile(outRareFile, outRare); outputNames.push_back(outRareFile); outputTypes["rarefy"].push_back(outRareFile); }
295 if (collect) { m->openOutputFile(outCollectFile, outCollect); outputNames.push_back(outCollectFile); outputTypes["phylodiv"].push_back(outCollectFile); }
297 int numLeafNodes = trees[i]->getNumLeaves();
299 //create a vector containing indexes of leaf nodes, randomize it, select nodes to send to calculator
300 vector<int> randomLeaf;
301 for (int j = 0; j < numLeafNodes; j++) {
302 if (m->inUsersGroups(trees[i]->tree[j].getGroup(), m->Groups) == true) { //is this a node from the group the user selected.
303 randomLeaf.push_back(j);
307 numLeafNodes = randomLeaf.size(); //reset the number of leaf nodes you are using
309 //each group, each sampling, if no rarefy iters = 1;
310 map<string, vector<float> > diversity;
312 //each group, each sampling, if no rarefy iters = 1;
313 map<string, vector<float> > sumDiversity;
315 //find largest group total
316 int largestGroup = 0;
317 for (int j = 0; j < m->Groups.size(); j++) {
318 if (tmap->seqsPerGroup[m->Groups[j]] > largestGroup) { largestGroup = tmap->seqsPerGroup[m->Groups[j]]; }
320 //initialize diversity
321 diversity[m->Groups[j]].resize(tmap->seqsPerGroup[m->Groups[j]]+1, 0.0); //numSampled
324 //initialize sumDiversity
325 sumDiversity[m->Groups[j]].resize(tmap->seqsPerGroup[m->Groups[j]]+1, 0.0);
328 //convert freq percentage to number
330 if (freq < 1.0) { increment = largestGroup * freq;
331 }else { increment = freq; }
333 //initialize sampling spots
334 set<int> numSampledList;
335 for(int k = 1; k <= largestGroup; k++){ if((k == 1) || (k % increment == 0)){ numSampledList.insert(k); } }
336 if(largestGroup % increment != 0){ numSampledList.insert(largestGroup); }
338 //add other groups ending points
339 for (int j = 0; j < m->Groups.size(); j++) {
340 if (numSampledList.count(diversity[m->Groups[j]].size()-1) == 0) { numSampledList.insert(diversity[m->Groups[j]].size()-1); }
343 #if defined (__APPLE__) || (__MACH__) || (linux) || (__linux)
345 driver(trees[i], diversity, sumDiversity, iters, increment, randomLeaf, numSampledList, outCollect, outSum, true);
348 vector<int> procIters;
350 int numItersPerProcessor = iters / processors;
352 //divide iters between processes
353 for (int h = 0; h < processors; h++) {
354 if(h == processors - 1){
355 numItersPerProcessor = iters - h * numItersPerProcessor;
357 procIters.push_back(numItersPerProcessor);
360 createProcesses(procIters, trees[i], diversity, sumDiversity, iters, increment, randomLeaf, numSampledList, outCollect, outSum);
362 }else{ //no need to paralellize if you dont want to rarefy
363 driver(trees[i], diversity, sumDiversity, iters, increment, randomLeaf, numSampledList, outCollect, outSum, true);
368 driver(trees[i], diversity, sumDiversity, iters, increment, randomLeaf, numSampledList, outCollect, outSum, true);
371 if (rarefy) { printData(numSampledList, sumDiversity, outRare, iters); }
375 if (m->control_pressed) { for (int i = 0; i < outputNames.size(); i++) { remove(outputNames[i].c_str()); } return 0; }
377 m->mothurOutEndLine();
378 m->mothurOut("Output File Names: "); m->mothurOutEndLine();
379 for (int i = 0; i < outputNames.size(); i++) { m->mothurOut(outputNames[i]); m->mothurOutEndLine(); }
380 m->mothurOutEndLine();
385 catch(exception& e) {
386 m->errorOut(e, "PhyloDiversityCommand", "execute");
390 //**********************************************************************************************************************
391 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){
393 #if defined (__APPLE__) || (__MACH__) || (linux) || (__linux)
396 vector<int> processIDS;
397 map< string, vector<float> >::iterator itSum;
399 //loop through and create all the processes you want
400 while (process != processors) {
404 processIDS.push_back(pid); //create map from line number to pid so you can append files in correct order later
407 driver(t, div, sumDiv, procIters[process], increment, randomLeaf, numSampledList, outCollect, outSum, false);
409 string outTemp = outputDir + toString(getpid()) + ".sumDiv.temp";
411 m->openOutputFile(outTemp, out);
413 //output the sumDIversity
414 for (itSum = sumDiv.begin(); itSum != sumDiv.end(); itSum++) {
415 out << itSum->first << '\t' << (itSum->second).size() << '\t';
416 for (int k = 0; k < (itSum->second).size(); k++) {
417 out << (itSum->second)[k] << '\t';
426 m->mothurOut("[ERROR]: unable to spawn the necessary processes."); m->mothurOutEndLine();
427 for (int i = 0; i < processIDS.size(); i++) { kill (processIDS[i], SIGINT); }
432 driver(t, div, sumDiv, procIters[0], increment, randomLeaf, numSampledList, outCollect, outSum, true);
434 //force parent to wait until all the processes are done
435 for (int i=0;i<(processors-1);i++) {
436 int temp = processIDS[i];
440 //get data created by processes
441 for (int i=0;i<(processors-1);i++) {
443 //input the sumDIversity
444 string inTemp = outputDir + toString(processIDS[i]) + ".sumDiv.temp";
446 m->openInputFile(inTemp, in);
448 //output the sumDIversity
449 for (int j = 0; j < sumDiv.size(); j++) {
453 in >> group >> size; m->gobble(in);
455 for (int k = 0; k < size; k++) {
459 sumDiv[group][k] += tempVal;
465 remove(inTemp.c_str());
473 catch(exception& e) {
474 m->errorOut(e, "PhyloDiversityCommand", "createProcesses");
478 //**********************************************************************************************************************
479 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){
481 int numLeafNodes = randomLeaf.size();
483 for (int l = 0; l < numIters; l++) {
484 random_shuffle(randomLeaf.begin(), randomLeaf.end());
487 map<string, int> counts;
488 map< string, set<int> > countedBranch;
489 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
491 for(int k = 0; k < numLeafNodes; k++){
493 if (m->control_pressed) { return 0; }
495 //calc branch length of randomLeaf k
496 vector<float> br = calcBranchLength(t, randomLeaf[k], countedBranch);
498 //for each group in the groups update the total branch length accounting for the names file
499 vector<string> groups = t->tree[randomLeaf[k]].getGroup();
501 for (int j = 0; j < groups.size(); j++) {
502 int numSeqsInGroupJ = 0;
503 map<string, int>::iterator it;
504 it = t->tree[randomLeaf[k]].pcount.find(groups[j]);
505 if (it != t->tree[randomLeaf[k]].pcount.end()) { //this leaf node contains seqs from group j
506 numSeqsInGroupJ = it->second;
509 if (numSeqsInGroupJ != 0) { div[groups[j]][(counts[groups[j]]+1)] = div[groups[j]][counts[groups[j]]] + br[j]; }
511 for (int s = (counts[groups[j]]+2); s <= (counts[groups[j]]+numSeqsInGroupJ); s++) {
512 div[groups[j]][s] = div[groups[j]][s-1]; //update counts, but don't add in redundant branch lengths
514 counts[groups[j]] += numSeqsInGroupJ;
519 //add this diversity to the sum
520 for (int j = 0; j < m->Groups.size(); j++) {
521 for (int g = 0; g < div[m->Groups[j]].size(); g++) {
522 sumDiv[m->Groups[j]][g] += div[m->Groups[j]][g];
527 if ((collect) && (l == 0) && doSumCollect) { printData(numSampledList, div, outCollect, 1); }
528 if ((summary) && (l == 0) && doSumCollect) { printSumData(div, outSum, 1); }
534 catch(exception& e) {
535 m->errorOut(e, "PhyloDiversityCommand", "driver");
540 //**********************************************************************************************************************
542 void PhyloDiversityCommand::printSumData(map< string, vector<float> >& div, ofstream& out, int numIters){
545 out << "Groups\tnumSampled\tphyloDiversity" << endl;
547 out.setf(ios::fixed, ios::floatfield); out.setf(ios::showpoint);
549 for (int j = 0; j < m->Groups.size(); j++) {
550 int numSampled = (div[m->Groups[j]].size()-1);
551 out << m->Groups[j] << '\t' << numSampled << '\t';
555 if (scale) { score = (div[m->Groups[j]][numSampled] / (float)numIters) / (float)numSampled; }
556 else { score = div[m->Groups[j]][numSampled] / (float)numIters; }
558 out << setprecision(4) << score << endl;
564 catch(exception& e) {
565 m->errorOut(e, "PhyloDiversityCommand", "printSumData");
569 //**********************************************************************************************************************
571 void PhyloDiversityCommand::printData(set<int>& num, map< string, vector<float> >& div, ofstream& out, int numIters){
574 out << "numSampled\t";
575 for (int i = 0; i < m->Groups.size(); i++) { out << m->Groups[i] << '\t'; }
578 out.setf(ios::fixed, ios::floatfield); out.setf(ios::showpoint);
580 for (set<int>::iterator it = num.begin(); it != num.end(); it++) {
581 int numSampled = *it;
583 out << numSampled << '\t';
585 for (int j = 0; j < m->Groups.size(); j++) {
586 if (numSampled < div[m->Groups[j]].size()) {
588 if (scale) { score = (div[m->Groups[j]][numSampled] / (float)numIters) / (float)numSampled; }
589 else { score = div[m->Groups[j]][numSampled] / (float)numIters; }
591 out << setprecision(4) << score << '\t';
592 }else { out << "NA" << '\t'; }
600 catch(exception& e) {
601 m->errorOut(e, "PhyloDiversityCommand", "printData");
605 //**********************************************************************************************************************
606 //need a vector of floats one branch length for every group the node represents.
607 vector<float> PhyloDiversityCommand::calcBranchLength(Tree* t, int leaf, map< string, set<int> >& counted){
610 //calc the branch length
611 //while you aren't at root
615 vector<string> groups = t->tree[leaf].getGroup();
616 sums.resize(groups.size(), 0.0);
618 map<string, map<int, double> > tempTotals; //maps node to total Branch Length
619 map< string, set<int> > tempCounted;
620 set<int>::iterator it;
623 if(t->tree[index].getBranchLength() != -1){
624 for (int k = 0; k < groups.size(); k++) {
625 sums[k] += abs(t->tree[index].getBranchLength());
626 counted[groups[k]].insert(index);
630 for (int k = 0; k < groups.size(); k++) {
631 tempTotals[groups[k]][index] = 0.0;
634 index = t->tree[index].getParent();
636 //while you aren't at root
637 while(t->tree[index].getParent() != -1){
639 if (m->control_pressed) { return sums; }
642 for (int k = 0; k < groups.size(); k++) {
643 map<string, int>::iterator itGroup = t->tree[index].pcount.find(groups[k]);
644 if (itGroup != t->tree[index].pcount.end()) { pcountSize++; }
646 //do both your chidren have have descendants from the users groups?
647 int lc = t->tree[index].getLChild();
648 int rc = t->tree[index].getRChild();
651 itGroup = t->tree[lc].pcount.find(groups[k]);
652 if (itGroup != t->tree[lc].pcount.end()) { LpcountSize++; }
655 itGroup = t->tree[rc].pcount.find(groups[k]);
656 if (itGroup != t->tree[rc].pcount.end()) { RpcountSize++; }
658 //if yes, add your childrens tempTotals
659 if ((LpcountSize != 0) && (RpcountSize != 0)) {
660 sums[k] += tempTotals[groups[k]][lc] + tempTotals[groups[k]][rc];
662 for (it = tempCounted[groups[k]].begin(); it != tempCounted[groups[k]].end(); it++) { counted[groups[k]].insert(*it); }
664 //cout << "added to total " << tempTotals[lc] << '\t' << tempTotals[rc] << endl;
665 if (t->tree[index].getBranchLength() != -1) {
666 if (counted[groups[k]].count(index) == 0) {
667 tempTotals[groups[k]][index] = abs(t->tree[index].getBranchLength());
668 tempCounted[groups[k]].insert(index);
670 tempTotals[groups[k]][index] = 0.0;
673 tempTotals[groups[k]][index] = 0.0;
675 }else { //if no, your tempTotal is your childrens temp totals + your branch length
676 tempTotals[groups[k]][index] = tempTotals[groups[k]][lc] + tempTotals[groups[k]][rc];
678 if (counted[groups[k]].count(index) == 0) {
679 tempTotals[groups[k]][index] += abs(t->tree[index].getBranchLength());
680 tempCounted[groups[k]].insert(index);
684 //cout << "temptotal = "<< tempTotals[i] << endl;
687 index = t->tree[index].getParent();
693 catch(exception& e) {
694 m->errorOut(e, "PhyloDiversityCommand", "calcBranchLength");
698 /*****************************************************************/
699 int PhyloDiversityCommand::readNamesFile() {
702 numUniquesInName = 0;
705 m->openInputFile(namefile, in);
707 string first, second;
708 map<string, string>::iterator itNames;
711 in >> first >> second; m->gobble(in);
715 itNames = m->names.find(first);
716 if (itNames == m->names.end()) {
717 m->names[first] = second;
719 //we need a list of names in your namefile to use above when removing extra seqs above so we don't remove them
720 vector<string> dupNames;
721 m->splitAtComma(second, dupNames);
723 for (int i = 0; i < dupNames.size(); i++) {
724 nameMap[dupNames[i]] = dupNames[i];
725 if ((groupfile == "") && (i != 0)) { tmap->addSeq(dupNames[i], "Group1"); }
727 }else { m->mothurOut(first + " has already been seen in namefile, disregarding names file."); m->mothurOutEndLine(); in.close(); m->names.clear(); namefile = ""; return 1; }
733 catch(exception& e) {
734 m->errorOut(e, "PhyloDiversityCommand", "readNamesFile");
739 //**********************************************************************************************************************