+ return 0;
+
+ }
+ catch(exception& e) {
+ m->errorOut(e, "PhyloDiversityCommand", "createProcesses");
+ exit(1);
+ }
+}
+//**********************************************************************************************************************
+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){
+ try {
+ int numLeafNodes = randomLeaf.size();
+
+ for (int l = 0; l < numIters; l++) {
+ random_shuffle(randomLeaf.begin(), randomLeaf.end());
+
+ //initialize counts
+ map<string, int> counts;
+ map< string, set<int> > countedBranch;
+ for (int j = 0; j < globaldata->Groups.size(); j++) { counts[globaldata->Groups[j]] = 0; countedBranch[globaldata->Groups[j]].insert(-2); } //add dummy index to initialize countedBranch sets
+
+ for(int k = 0; k < numLeafNodes; k++){
+
+ if (m->control_pressed) { return 0; }
+
+ //calc branch length of randomLeaf k
+ vector<float> br = calcBranchLength(t, randomLeaf[k], countedBranch);
+
+ //for each group in the groups update the total branch length accounting for the names file
+ vector<string> groups = t->tree[randomLeaf[k]].getGroup();
+
+ for (int j = 0; j < groups.size(); j++) {
+ int numSeqsInGroupJ = 0;
+ map<string, int>::iterator it;
+ it = t->tree[randomLeaf[k]].pcount.find(groups[j]);
+ if (it != t->tree[randomLeaf[k]].pcount.end()) { //this leaf node contains seqs from group j
+ numSeqsInGroupJ = it->second;
+ }
+
+ if (numSeqsInGroupJ != 0) { div[groups[j]][(counts[groups[j]]+1)] = div[groups[j]][counts[groups[j]]] + br[j]; }
+
+ for (int s = (counts[groups[j]]+2); s <= (counts[groups[j]]+numSeqsInGroupJ); s++) {
+ div[groups[j]][s] = div[groups[j]][s-1]; //update counts, but don't add in redundant branch lengths
+ }
+ counts[groups[j]] += numSeqsInGroupJ;
+ }
+ }
+
+ if (rarefy) {
+ //add this diversity to the sum
+ for (int j = 0; j < globaldata->Groups.size(); j++) {
+ for (int g = 0; g < div[globaldata->Groups[j]].size(); g++) {
+ sumDiv[globaldata->Groups[j]][g] += div[globaldata->Groups[j]][g];
+ }
+ }
+ }
+
+ if ((collect) && (l == 0) && doSumCollect) { printData(numSampledList, div, outCollect, 1); }
+ if ((summary) && (l == 0) && doSumCollect) { printSumData(div, outSum, 1); }
+ }
+
+ return 0;
+
+ }
+ catch(exception& e) {
+ m->errorOut(e, "PhyloDiversityCommand", "driver");
+ exit(1);
+ }
+}
+
+//**********************************************************************************************************************
+
+void PhyloDiversityCommand::printSumData(map< string, vector<float> >& div, ofstream& out, int numIters){
+ try {
+
+ out << "Groups\tnumSampled\tphyloDiversity" << endl;
+
+ out.setf(ios::fixed, ios::floatfield); out.setf(ios::showpoint);
+
+ for (int j = 0; j < globaldata->Groups.size(); j++) {
+ int numSampled = (div[globaldata->Groups[j]].size()-1);
+ out << globaldata->Groups[j] << '\t' << numSampled << '\t';
+
+
+ float score;
+ if (scale) { score = (div[globaldata->Groups[j]][numSampled] / (float)numIters) / (float)numSampled; }
+ else { score = div[globaldata->Groups[j]][numSampled] / (float)numIters; }
+
+ out << setprecision(4) << score << endl;
+ }
+
+ out.close();
+
+ }
+ catch(exception& e) {
+ m->errorOut(e, "PhyloDiversityCommand", "printSumData");
+ exit(1);
+ }
+}
+//**********************************************************************************************************************
+
+void PhyloDiversityCommand::printData(set<int>& num, map< string, vector<float> >& div, ofstream& out, int numIters){