]> git.donarmstrong.com Git - mothur.git/blob - phylodiversitycommand.cpp
added get.otulabels and remove.otulabels commands per user request
[mothur.git] / phylodiversitycommand.cpp
1 /*
2  *  phylodiversitycommand.cpp
3  *  Mothur
4  *
5  *  Created by westcott on 4/30/10.
6  *  Copyright 2010 Schloss Lab. All rights reserved.
7  *
8  */
9
10 #include "phylodiversitycommand.h"
11 #include "treereader.h"
12
13 //**********************************************************************************************************************
14 vector<string> PhyloDiversityCommand::setParameters(){  
15         try {
16
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);
30                 
31                 vector<string> myArray;
32                 for (int i = 0; i < parameters.size(); i++) {   myArray.push_back(parameters[i].name);          }
33                 return myArray;
34         }
35         catch(exception& e) {
36                 m->errorOut(e, "PhyloDiversityCommand", "setParameters");
37                 exit(1);
38         }
39 }
40 //**********************************************************************************************************************
41 string PhyloDiversityCommand::getHelpString(){  
42         try {
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";
57                 return helpString;
58         }
59         catch(exception& e) {
60                 m->errorOut(e, "PhyloDiversityCommand", "getHelpString");
61                 exit(1);
62         }
63 }
64
65
66 //**********************************************************************************************************************
67 PhyloDiversityCommand::PhyloDiversityCommand(){ 
68         try {
69                 abort = true; calledHelp = true; 
70                 setParameters();
71                 vector<string> tempOutNames;
72                 outputTypes["phylodiv"] = tempOutNames;
73                 outputTypes["rarefy"] = tempOutNames;
74                 outputTypes["summary"] = tempOutNames;
75         }
76         catch(exception& e) {
77                 m->errorOut(e, "PhyloDiversityCommand", "PhyloDiversityCommand");
78                 exit(1);
79         }
80 }
81 //**********************************************************************************************************************
82 PhyloDiversityCommand::PhyloDiversityCommand(string option)  {
83         try {
84                 abort = false; calledHelp = false;   
85                 
86                 //allow user to run help
87                 if(option == "help") { help(); abort = true; calledHelp = true; }
88                 else if(option == "citation") { citation(); abort = true; calledHelp = true;}
89                 
90                 else {
91                         vector<string> myArray = setParameters();;
92                         
93                         OptionParser parser(option);
94                         map<string,string> parameters = parser.getParameters();
95                         map<string,string>::iterator it;
96                         
97                         ValidParameters validParameter;
98                 
99                         //check to make sure all parameters are valid for command
100                         for (map<string,string>::iterator it = parameters.begin(); it != parameters.end(); it++) { 
101                                 if (validParameter.isValidParameter(it->first, myArray, it->second) != true) {  abort = true;  }
102                         }
103                         
104                         //initialize outputTypes
105                         vector<string> tempOutNames;
106                         outputTypes["phylodiv"] = tempOutNames;
107                         outputTypes["rarefy"] = tempOutNames;
108                         outputTypes["summary"] = tempOutNames;
109                         
110                         //if the user changes the input directory command factory will send this info to us in the output parameter 
111                         string inputDir = validParameter.validFile(parameters, "inputdir", false);              
112                         if (inputDir == "not found"){   inputDir = "";          }
113                         else {
114                                 string path;
115                                 it = parameters.find("tree");
116                                 //user has given a template file
117                                 if(it != parameters.end()){ 
118                                         path = m->hasPath(it->second);
119                                         //if the user has not given a path then, add inputdir. else leave path alone.
120                                         if (path == "") {       parameters["tree"] = inputDir + it->second;             }
121                                 }
122                                 
123                                 it = parameters.find("group");
124                                 //user has given a template file
125                                 if(it != parameters.end()){ 
126                                         path = m->hasPath(it->second);
127                                         //if the user has not given a path then, add inputdir. else leave path alone.
128                                         if (path == "") {       parameters["group"] = inputDir + it->second;            }
129                                 }
130                                 
131                                 it = parameters.find("name");
132                                 //user has given a template file
133                                 if(it != parameters.end()){ 
134                                         path = m->hasPath(it->second);
135                                         //if the user has not given a path then, add inputdir. else leave path alone.
136                                         if (path == "") {       parameters["name"] = inputDir + it->second;             }
137                                 }
138                         }
139                         
140                         //check for required parameters
141                         treefile = validParameter.validFile(parameters, "tree", true);
142                         if (treefile == "not open") { treefile = ""; abort = true; }
143                         else if (treefile == "not found") {                             
144                                 //if there is a current design file, use it
145                                 treefile = m->getTreeFile(); 
146                                 if (treefile != "") { m->mothurOut("Using " + treefile + " as input file for the tree parameter."); m->mothurOutEndLine(); }
147                                 else {  m->mothurOut("You have no current tree file and the tree parameter is required."); m->mothurOutEndLine(); abort = true; }                                                               
148                         }else { m->setTreeFile(treefile); }     
149                         
150                         //check for required parameters
151                         groupfile = validParameter.validFile(parameters, "group", true);
152                         if (groupfile == "not open") { groupfile = ""; abort = true; }
153                         else if (groupfile == "not found") { groupfile = ""; }
154                         else { m->setGroupFile(groupfile); }
155                         
156                         namefile = validParameter.validFile(parameters, "name", true);
157                         if (namefile == "not open") { namefile = ""; abort = true; }
158                         else if (namefile == "not found") { namefile = ""; }
159                         else { m->setNameFile(namefile); }
160                         
161                         outputDir = validParameter.validFile(parameters, "outputdir", false);           if (outputDir == "not found"){  outputDir = m->hasPath(treefile);       }
162                         
163                         string temp;
164                         temp = validParameter.validFile(parameters, "freq", false);                     if (temp == "not found") { temp = "100"; }
165                         m->mothurConvert(temp, freq); 
166                         
167                         temp = validParameter.validFile(parameters, "iters", false);                    if (temp == "not found") { temp = "1000"; }
168                         m->mothurConvert(temp, iters); 
169                         
170                         temp = validParameter.validFile(parameters, "rarefy", false);                   if (temp == "not found") { temp = "F"; }
171                         rarefy = m->isTrue(temp);
172                         if (!rarefy) { iters = 1;  }
173                         
174                         temp = validParameter.validFile(parameters, "summary", false);                  if (temp == "not found") { temp = "T"; }
175                         summary = m->isTrue(temp);
176                         
177                         temp = validParameter.validFile(parameters, "scale", false);                    if (temp == "not found") { temp = "F"; }
178                         scale = m->isTrue(temp);
179                         
180                         temp = validParameter.validFile(parameters, "collect", false);                  if (temp == "not found") { temp = "F"; }
181                         collect = m->isTrue(temp);
182                         
183                         temp = validParameter.validFile(parameters, "processors", false);       if (temp == "not found"){       temp = m->getProcessors();      }
184                         m->setProcessors(temp);
185                         m->mothurConvert(temp, processors); 
186                         
187                         groups = validParameter.validFile(parameters, "groups", false);                 
188                         if (groups == "not found") { groups = "";  }
189                         else { 
190                                 m->splitAtDash(groups, Groups);
191                                 m->setGroups(Groups);
192                         }
193                         
194                         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; }
195                         
196                         if (namefile == "") {
197                                 vector<string> files; files.push_back(treefile);
198                                 parser.getNameFile(files);
199                         }
200                 }
201                 
202         }
203         catch(exception& e) {
204                 m->errorOut(e, "PhyloDiversityCommand", "PhyloDiversityCommand");
205                 exit(1);
206         }                       
207 }
208 //**********************************************************************************************************************
209
210 int PhyloDiversityCommand::execute(){
211         try {
212                 
213                 if (abort == true) { if (calledHelp) { return 0; }  return 2;   }
214                 
215         int start = time(NULL);
216         
217                 m->setTreeFile(treefile);
218         TreeReader* reader = new TreeReader(treefile, groupfile, namefile);
219         vector<Tree*> trees = reader->getTrees();
220         tmap = trees[0]->getTreeMap();
221         delete reader;
222
223                 SharedUtil util;
224                 vector<string> mGroups = m->getGroups();
225                 vector<string> tGroups = tmap->getNamesOfGroups();
226                 util.setGroups(mGroups, tGroups, "phylo.diversity");    //sets the groups the user wants to analyze
227                 
228                 //incase the user had some mismatches between the tree and group files we don't want group xxx to be analyzed
229                 for (int i = 0; i < mGroups.size(); i++) { if (mGroups[i] == "xxx") { mGroups.erase(mGroups.begin()+i);  break; }  }
230                 m->setGroups(mGroups);
231                  
232                 vector<string> outputNames;
233                 
234                 //for each of the users trees
235                 for(int i = 0; i < trees.size(); i++) {
236                 
237                         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; }
238                         
239                         ofstream outSum, outRare, outCollect;
240                         string outSumFile = outputDir + m->getRootName(m->getSimpleName(treefile))  + toString(i+1) + ".phylodiv.summary";
241                         string outRareFile = outputDir + m->getRootName(m->getSimpleName(treefile))  + toString(i+1) + ".phylodiv.rarefaction";
242                         string outCollectFile = outputDir + m->getRootName(m->getSimpleName(treefile))  + toString(i+1) + ".phylodiv";
243                         
244                         if (summary)    { m->openOutputFile(outSumFile, outSum); outputNames.push_back(outSumFile);             outputTypes["summary"].push_back(outSumFile);                   }
245                         if (rarefy)             { m->openOutputFile(outRareFile, outRare); outputNames.push_back(outRareFile);  outputTypes["rarefy"].push_back(outRareFile);                   }
246                         if (collect)    { m->openOutputFile(outCollectFile, outCollect); outputNames.push_back(outCollectFile);  outputTypes["phylodiv"].push_back(outCollectFile);  }
247                         
248                         int numLeafNodes = trees[i]->getNumLeaves();
249                         
250                         //create a vector containing indexes of leaf nodes, randomize it, select nodes to send to calculator
251                         vector<int> randomLeaf;
252                         for (int j = 0; j < numLeafNodes; j++) {  
253                                 if (m->inUsersGroups(trees[i]->tree[j].getGroup(), mGroups) == true) { //is this a node from the group the user selected.
254                                         randomLeaf.push_back(j); 
255                                 }
256                         }
257                         
258                         numLeafNodes = randomLeaf.size();  //reset the number of leaf nodes you are using 
259                         
260                         //each group, each sampling, if no rarefy iters = 1;
261                         map<string, vector<float> > diversity;
262                         
263                         //each group, each sampling, if no rarefy iters = 1;
264                         map<string, vector<float> > sumDiversity;
265                         
266                         //find largest group total 
267                         int largestGroup = 0;
268                         for (int j = 0; j < mGroups.size(); j++) {  
269                                 if (tmap->seqsPerGroup[mGroups[j]] > largestGroup) { largestGroup = tmap->seqsPerGroup[mGroups[j]]; }
270                                 
271                                 //initialize diversity
272                                 diversity[mGroups[j]].resize(tmap->seqsPerGroup[mGroups[j]]+1, 0.0);            //numSampled
273                                                                                                                                                                                                                         //groupA                0.0                     0.0
274                                                                                                                                                                                                                         
275                                 //initialize sumDiversity
276                                 sumDiversity[mGroups[j]].resize(tmap->seqsPerGroup[mGroups[j]]+1, 0.0);
277                         }       
278
279                         //convert freq percentage to number
280                         int increment = 100;
281                         if (freq < 1.0) {  increment = largestGroup * freq;  
282                         }else { increment = freq;  }
283                         
284                         //initialize sampling spots
285                         set<int> numSampledList;
286                         for(int k = 1; k <= largestGroup; k++){  if((k == 1) || (k % increment == 0)){  numSampledList.insert(k); }   }
287                         if(largestGroup % increment != 0){      numSampledList.insert(largestGroup);   }
288                         
289                         //add other groups ending points
290                         for (int j = 0; j < mGroups.size(); j++) {  
291                                 if (numSampledList.count(diversity[mGroups[j]].size()-1) == 0) {  numSampledList.insert(diversity[mGroups[j]].size()-1); }
292                         }
293                         
294                         #if defined (__APPLE__) || (__MACH__) || (linux) || (__linux) || (__linux__) || (__unix__) || (__unix)
295                                 if(processors == 1){
296                                         driver(trees[i], diversity, sumDiversity, iters, increment, randomLeaf, numSampledList, outCollect, outSum, true);      
297                                 }else{
298                                         if (rarefy) {
299                                                 vector<int> procIters;
300                                                 
301                                                 int numItersPerProcessor = iters / processors;
302                                                 
303                                                 //divide iters between processes
304                                                 for (int h = 0; h < processors; h++) {
305                                                         if(h == processors - 1){
306                                                                 numItersPerProcessor = iters - h * numItersPerProcessor;
307                                                         }
308                                                         procIters.push_back(numItersPerProcessor);
309                                                 }
310                                                 
311                                                 createProcesses(procIters, trees[i], diversity, sumDiversity, iters, increment, randomLeaf, numSampledList, outCollect, outSum); 
312                                                 
313                                         }else{ //no need to paralellize if you dont want to rarefy
314                                                 driver(trees[i], diversity, sumDiversity, iters, increment, randomLeaf, numSampledList, outCollect, outSum, true);      
315                                         }
316                                 }
317
318                         #else
319                                 driver(trees[i], diversity, sumDiversity, iters, increment, randomLeaf, numSampledList, outCollect, outSum, true);      
320                         #endif
321
322                         if (rarefy) {   printData(numSampledList, sumDiversity, outRare, iters);        }
323                 }
324                 
325         
326                 if (m->control_pressed) { for (int i = 0; i < outputNames.size(); i++) {        m->mothurRemove(outputNames[i]);        } return 0; }
327
328         m->mothurOut("It took " + toString(time(NULL) - start) + " secs to run unifrac.unweighted."); m->mothurOutEndLine();
329
330         
331                 m->mothurOutEndLine();
332                 m->mothurOut("Output File Names: "); m->mothurOutEndLine();
333                 for (int i = 0; i < outputNames.size(); i++) {  m->mothurOut(outputNames[i]); m->mothurOutEndLine();    }
334                 m->mothurOutEndLine();
335
336                 
337                 return 0;
338         }
339         catch(exception& e) {
340                 m->errorOut(e, "PhyloDiversityCommand", "execute");
341                 exit(1);
342         }
343 }
344 //**********************************************************************************************************************
345 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){
346         try {
347                 #if defined (__APPLE__) || (__MACH__) || (linux) || (__linux) || (__linux__) || (__unix__) || (__unix)
348                 int process = 1;
349                 
350                 vector<int> processIDS;
351                 map< string, vector<float> >::iterator itSum;
352                 
353                 //loop through and create all the processes you want
354                 while (process != processors) {
355                         int pid = fork();
356                         
357                         if (pid > 0) {
358                                 processIDS.push_back(pid);  //create map from line number to pid so you can append files in correct order later
359                                 process++;
360                         }else if (pid == 0){
361                                 driver(t, div, sumDiv, procIters[process], increment, randomLeaf, numSampledList, outCollect, outSum, false);
362                                 
363                                 string outTemp = outputDir + toString(getpid()) + ".sumDiv.temp";
364                                 ofstream out;
365                                 m->openOutputFile(outTemp, out);
366                                 
367                                 //output the sumDIversity
368                                 for (itSum = sumDiv.begin(); itSum != sumDiv.end(); itSum++) {
369                                         out << itSum->first << '\t' << (itSum->second).size() << '\t';
370                                         for (int k = 0; k < (itSum->second).size(); k++) { 
371                                                 out << (itSum->second)[k] << '\t';
372                                         }
373                                         out << endl;
374                                 }
375                                 
376                                 out.close();
377                                 
378                                 exit(0);
379                         }else { 
380                                 m->mothurOut("[ERROR]: unable to spawn the necessary processes."); m->mothurOutEndLine(); 
381                                 for (int i = 0; i < processIDS.size(); i++) { kill (processIDS[i], SIGINT); }
382                                 exit(0);
383                         }
384                 }
385                 
386                 driver(t, div, sumDiv, procIters[0], increment, randomLeaf, numSampledList, outCollect, outSum, true);
387                 
388                 //force parent to wait until all the processes are done
389                 for (int i=0;i<(processors-1);i++) { 
390                         int temp = processIDS[i];
391                         wait(&temp);
392                 }
393                 
394                 //get data created by processes
395                 for (int i=0;i<(processors-1);i++) { 
396                         
397                         //input the sumDIversity
398                         string inTemp = outputDir + toString(processIDS[i]) + ".sumDiv.temp";
399                         ifstream in;
400                         m->openInputFile(inTemp, in);
401                                 
402                         //output the sumDIversity
403                         for (int j = 0; j < sumDiv.size(); j++) { 
404                                 string group = "";
405                                 int size = 0;
406                                 
407                                 in >> group >> size; m->gobble(in);
408                                 
409                                 for (int k = 0; k < size; k++) { 
410                                         float tempVal;
411                                         in >> tempVal;
412                                         
413                                         sumDiv[group][k] += tempVal;
414                                 }
415                                 m->gobble(in);
416                         }
417                                 
418                         in.close();
419                         m->mothurRemove(inTemp);
420                 }
421                 
422 #endif
423
424         return 0;               
425         
426         }
427         catch(exception& e) {
428                 m->errorOut(e, "PhyloDiversityCommand", "createProcesses");
429                 exit(1);
430         }
431 }
432 //**********************************************************************************************************************
433 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){
434         try {
435                 int numLeafNodes = randomLeaf.size();
436                 vector<string> mGroups = m->getGroups();
437         
438         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.
439         
440                 for (int l = 0; l < numIters; l++) {
441                                 random_shuffle(randomLeaf.begin(), randomLeaf.end());
442             cout << l << endl;
443                                 //initialize counts
444                                 map<string, int> counts;
445                 vector< map<string, bool> > countedBranch;
446                 for (int i = 0; i < t->getNumNodes(); i++) {
447                     map<string, bool> temp;
448                     for (int j = 0; j < mGroups.size(); j++) { temp[mGroups[j]] = false; }
449                     countedBranch.push_back(temp);
450                 }
451             
452                                 for (int j = 0; j < mGroups.size(); j++) {  counts[mGroups[j]] = 0;   }  
453                                 
454                                 for(int k = 0; k < numLeafNodes; k++){
455                                                 
456                                         if (m->control_pressed) { return 0; }
457                                         
458                                         //calc branch length of randomLeaf k
459                     vector<float> br = calcBranchLength(t, randomLeaf[k], countedBranch, rootForGroup);
460                         
461                                         //for each group in the groups update the total branch length accounting for the names file
462                                         vector<string> groups = t->tree[randomLeaf[k]].getGroup();
463                                         
464                                         for (int j = 0; j < groups.size(); j++) {
465                         
466                         if (m->inUsersGroups(groups[j], mGroups)) {
467                             int numSeqsInGroupJ = 0;
468                             map<string, int>::iterator it;
469                             it = t->tree[randomLeaf[k]].pcount.find(groups[j]);
470                             if (it != t->tree[randomLeaf[k]].pcount.end()) { //this leaf node contains seqs from group j
471                                 numSeqsInGroupJ = it->second;
472                             }
473                             
474                             if (numSeqsInGroupJ != 0) { div[groups[j]][(counts[groups[j]]+1)] = div[groups[j]][counts[groups[j]]] + br[j];  }
475                             
476                             for (int s = (counts[groups[j]]+2); s <= (counts[groups[j]]+numSeqsInGroupJ); s++) {
477                                 div[groups[j]][s] = div[groups[j]][s-1];  //update counts, but don't add in redundant branch lengths
478                             }
479                             counts[groups[j]] += numSeqsInGroupJ;
480                         }
481                                         }
482                                 }
483                                 
484                                 if (rarefy) {
485                                         //add this diversity to the sum
486                                         for (int j = 0; j < mGroups.size(); j++) {  
487                                                 for (int g = 0; g < div[mGroups[j]].size(); g++) {
488                                                         sumDiv[mGroups[j]][g] += div[mGroups[j]][g];
489                                                 }
490                                         }
491                                 }
492                                 
493                                 if ((collect) && (l == 0) && doSumCollect) {  printData(numSampledList, div, outCollect, 1);  }
494                                 if ((summary) && (l == 0) && doSumCollect) {  printSumData(div, outSum, 1);  }
495                         }
496                         
497                         return 0;
498
499         }
500         catch(exception& e) {
501                 m->errorOut(e, "PhyloDiversityCommand", "driver");
502                 exit(1);
503         }
504 }
505
506 //**********************************************************************************************************************
507
508 void PhyloDiversityCommand::printSumData(map< string, vector<float> >& div, ofstream& out, int numIters){
509         try {
510                 
511                 out << "Groups\tnumSampled\tphyloDiversity" << endl;
512                 
513                 out.setf(ios::fixed, ios::floatfield); out.setf(ios::showpoint);
514                 
515                 vector<string> mGroups = m->getGroups();
516                 for (int j = 0; j < mGroups.size(); j++) {
517                         int numSampled = (div[mGroups[j]].size()-1);
518                         out << mGroups[j] << '\t' << numSampled << '\t';
519                 
520                          
521                         float score;
522                         if (scale)      {  score = (div[mGroups[j]][numSampled] / (float)numIters) / (float)numSampled; }
523                         else            {       score = div[mGroups[j]][numSampled] / (float)numIters;  }
524                                 
525                         out << setprecision(4) << score << endl;
526                 }
527                                         
528                 out.close();
529                 
530         }
531         catch(exception& e) {
532                 m->errorOut(e, "PhyloDiversityCommand", "printSumData");
533                 exit(1);
534         }
535 }
536 //**********************************************************************************************************************
537
538 void PhyloDiversityCommand::printData(set<int>& num, map< string, vector<float> >& div, ofstream& out, int numIters){
539         try {
540                 
541                 out << "numSampled\t";
542                 vector<string> mGroups = m->getGroups();
543                 for (int i = 0; i < mGroups.size(); i++) { out << mGroups[i] << '\t';  }
544                 out << endl;
545                 
546                 out.setf(ios::fixed, ios::floatfield); out.setf(ios::showpoint);
547                 
548                 for (set<int>::iterator it = num.begin(); it != num.end(); it++) {  
549                         int numSampled = *it;
550                         
551                         out << numSampled << '\t';  
552                 
553                         for (int j = 0; j < mGroups.size(); j++) {
554                                 if (numSampled < div[mGroups[j]].size()) { 
555                                         float score;
556                                         if (scale)      {  score = (div[mGroups[j]][numSampled] / (float)numIters) / (float)numSampled; }
557                                         else            {       score = div[mGroups[j]][numSampled] / (float)numIters;  }
558
559                                         out << setprecision(4) << score << '\t';
560                                 }else { out << "NA" << '\t'; }
561                         }
562                         out << endl;
563                 }
564                 
565                 out.close();
566                 
567         }
568         catch(exception& e) {
569                 m->errorOut(e, "PhyloDiversityCommand", "printData");
570                 exit(1);
571         }
572 }
573 //**********************************************************************************************************************
574 //need a vector of floats one branch length for every group the node represents.
575 vector<float> PhyloDiversityCommand::calcBranchLength(Tree* t, int leaf, vector< map<string, bool> >& counted, map<string, int> roots){
576         try {
577         
578                 //calc the branch length
579                 //while you aren't at root
580                 vector<float> sums; 
581                 int index = leaf;
582                 
583                 vector<string> groups = t->tree[leaf].getGroup();
584                 sums.resize(groups.size(), 0.0);
585                 
586         
587         //you are a leaf
588                 if(t->tree[index].getBranchLength() != -1){     
589                         for (int k = 0; k < groups.size(); k++) { 
590                                 sums[k] += abs(t->tree[index].getBranchLength());       
591                         }
592                 }
593         
594         
595         index = t->tree[index].getParent();     
596         
597                 //while you aren't at root
598                 while(t->tree[index].getParent() != -1){
599             
600                         if (m->control_pressed) {  return sums; }
601                         
602                         for (int k = 0; k < groups.size(); k++) {
603                 
604                 if (index >= roots[groups[k]]) { counted[index][groups[k]] = true; } //if you are at this groups "root", then say we are done
605                 
606                 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
607                     if (t->tree[index].getBranchLength() != -1) {
608                         sums[k] += abs(t->tree[index].getBranchLength());
609                     }
610                     counted[index][groups[k]] = true;
611                 }
612             }
613             index = t->tree[index].getParent(); 
614         }
615         
616                 return sums;
617         
618         }
619         catch(exception& e) {
620                 m->errorOut(e, "PhyloDiversityCommand", "calcBranchLength");
621                 exit(1);
622         }
623 }
624 //**********************************************************************************************************************
625 map<string, int> PhyloDiversityCommand::getRootForGroups(Tree* t){
626         try {
627                 map<string, int> roots; //maps group to root for group, may not be root of tree
628                 map<string, bool> done;
629        
630                 //initialize root for all groups to -1
631                 for (int k = 0; k < (t->getTreeMap())->getNamesOfGroups().size(); k++) { done[(t->getTreeMap())->getNamesOfGroups()[k]] = false; }
632         
633         for (int i = 0; i < t->getNumLeaves(); i++) {
634             
635             vector<string> groups = t->tree[i].getGroup();
636             
637             int index = t->tree[i].getParent();
638             
639             for (int j = 0; j < groups.size(); j++) {
640                 
641                 if (done[groups[j]] == false) { //we haven't found the root for this group yet
642                     
643                     done[groups[j]] = true;
644                     roots[groups[j]] = i; //set root to self to start
645                     
646                     //while you aren't at root
647                     while(t->tree[index].getParent() != -1){
648                         
649                         if (m->control_pressed) {  return roots; }
650                         
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();
654                         
655                         int LpcountSize = 0;
656                         map<string, int>:: iterator itGroup = t->tree[lc].pcount.find(groups[j]);
657                         if (itGroup != t->tree[lc].pcount.end()) { LpcountSize++;  } 
658                         
659                         int RpcountSize = 0;
660                         itGroup = t->tree[rc].pcount.find(groups[j]);
661                         if (itGroup != t->tree[rc].pcount.end()) { RpcountSize++;  } 
662                         
663                         if ((LpcountSize != 0) && (RpcountSize != 0)) { //possible root
664                             roots[groups[j]] = index;
665                         }else { ;}
666                         
667                         index = t->tree[index].getParent();     
668                     }
669                 }
670             }
671         }
672         
673         return roots;
674         
675         }
676         catch(exception& e) {
677                 m->errorOut(e, "PhyloDiversityCommand", "getRootForGroups");
678                 exit(1);
679         }
680 }
681 //**********************************************************************************************************************
682
683
684