5 * Created by Sarah Westcott on 2/9/09.
6 * Copyright 2009 Schloss Lab UMASS Amherst. All rights reserved.
10 #include "unweighted.h"
12 /**************************************************************************************************/
14 EstOutput Unweighted::getValues(Tree* t, int p, string o) {
19 TreeMap* tmap = t->getTreeMap();
21 //if the users enters no groups then give them the score of all groups
22 int numGroups = m->getNumGroups();
24 //calculate number of comparsions
26 vector< vector<string> > namesOfGroupCombos;
27 for (int r=0; r<numGroups; r++) {
28 for (int l = 0; l < r; l++) {
30 vector<string> groups; groups.push_back((m->getGroups())[r]); groups.push_back((m->getGroups())[l]);
31 namesOfGroupCombos.push_back(groups);
36 vector<string> groups;
38 //get score for all users groups
39 for (int i = 0; i < (tmap->getNamesOfGroups()).size(); i++) {
40 if ((tmap->getNamesOfGroups())[i] != "xxx") {
41 groups.push_back((tmap->getNamesOfGroups())[i]);
44 namesOfGroupCombos.push_back(groups);
46 for (int i = 0; i < m->getNumGroups(); i++) {
47 groups.push_back((m->getGroups())[i]);
49 namesOfGroupCombos.push_back(groups);
53 #if defined (__APPLE__) || (__MACH__) || (linux) || (__linux) || (__linux__) || (__unix__) || (__unix)
55 data = driver(t, namesOfGroupCombos, 0, namesOfGroupCombos.size(), tmap);
57 int numPairs = namesOfGroupCombos.size();
59 int numPairsPerProcessor = numPairs / processors;
61 for (int i = 0; i < processors; i++) {
62 int startPos = i * numPairsPerProcessor;
64 if(i == processors - 1){
65 numPairsPerProcessor = numPairs - i * numPairsPerProcessor;
68 lines.push_back(linePair(startPos, numPairsPerProcessor));
70 data = createProcesses(t, namesOfGroupCombos, tmap);
74 data = driver(t, namesOfGroupCombos, 0, namesOfGroupCombos.size(), tmap);
80 m->errorOut(e, "Unweighted", "getValues");
84 /**************************************************************************************************/
86 EstOutput Unweighted::createProcesses(Tree* t, vector< vector<string> > namesOfGroupCombos, TreeMap* tmap) {
88 #if defined (__APPLE__) || (__MACH__) || (linux) || (__linux) || (__linux__) || (__unix__) || (__unix)
90 vector<int> processIDS;
94 //loop through and create all the processes you want
95 while (process != processors) {
99 processIDS.push_back(pid); //create map from line number to pid so you can append files in correct order later
103 myresults = driver(t, namesOfGroupCombos, lines[process].start, lines[process].num, tmap);
105 if (m->control_pressed) { exit(0); }
107 m->mothurOut("Merging results."); m->mothurOutEndLine();
109 //pass numSeqs to parent
111 string tempFile = outputDir + toString(getpid()) + ".unweighted.results.temp";
112 m->openOutputFile(tempFile, out);
113 out << myresults.size() << endl;
114 for (int i = 0; i < myresults.size(); i++) { out << myresults[i] << '\t'; } out << endl;
119 m->mothurOut("[ERROR]: unable to spawn the necessary processes."); m->mothurOutEndLine();
120 for (int i = 0; i < processIDS.size(); i++) { kill (processIDS[i], SIGINT); }
125 results = driver(t, namesOfGroupCombos, lines[0].start, lines[0].num, tmap);
127 //force parent to wait until all the processes are done
128 for (int i=0;i<(processors-1);i++) {
129 int temp = processIDS[i];
133 if (m->control_pressed) { return results; }
135 //get data created by processes
136 for (int i=0;i<(processors-1);i++) {
138 string s = outputDir + toString(processIDS[i]) + ".unweighted.results.temp";
139 m->openInputFile(s, in);
144 in >> num; m->gobble(in);
146 if (m->control_pressed) { break; }
149 for (int j = 0; j < num; j++) {
151 results.push_back(w);
159 m->mothurOut("DONE."); m->mothurOutEndLine(); m->mothurOutEndLine();
164 catch(exception& e) {
165 m->errorOut(e, "Unweighted", "createProcesses");
169 /**************************************************************************************************/
170 EstOutput Unweighted::driver(Tree* t, vector< vector<string> > namesOfGroupCombos, int start, int num, TreeMap* tmap) {
174 EstOutput results; results.resize(num);
178 int twentyPercent = (total * 0.20);
179 if (twentyPercent == 0) { twentyPercent = 1; }
182 for (int h = start; h < (start+num); h++) {
184 if (m->control_pressed) { return results; }
186 double UniqueBL=0.0000; //a branch length is unique if it's chidren are from the same group
187 double totalBL = 0.00; //all branch lengths
188 double UW = 0.00; //Unweighted Value = UniqueBL / totalBL;
190 //find a node that belongs to one of the groups in this combo
191 int nodeBelonging = -1;
192 for (int g = 0; g < namesOfGroupCombos[h].size(); g++) {
193 if (t->groupNodeInfo[namesOfGroupCombos[h][g]].size() != 0) { nodeBelonging = t->groupNodeInfo[namesOfGroupCombos[h][g]][0]; break; }
197 if (nodeBelonging == -1) {
198 m->mothurOut("[WARNING]: cannot find a nodes in the tree from grouping ");
199 for (int g = 0; g < namesOfGroupCombos[h].size()-1; g++) { m->mothurOut(namesOfGroupCombos[h][g] + "-"); }
200 m->mothurOut(namesOfGroupCombos[h][namesOfGroupCombos[h].size()-1]);
201 m->mothurOut(", skipping."); m->mothurOutEndLine(); results[count] = UW;
203 //cout << "trying to get root" << endl;
204 //if including the root this clears rootForGrouping[namesOfGroupCombos[h]]
205 getRoot(t, nodeBelonging, namesOfGroupCombos[h]);
206 //cout << "here" << endl;
207 for(int i=0;i<t->getNumNodes();i++){
209 if (m->control_pressed) { return data; }
211 //pcountSize = 0, they are from a branch that is entirely from a group the user doesn't want
212 //pcountSize = 2, not unique to one group
213 //pcountSize = 1, unique to one group
216 for (int j = 0; j < namesOfGroupCombos[h].size(); j++) {
217 map<string, int>::iterator itGroup = t->tree[i].pcount.find(namesOfGroupCombos[h][j]);
218 if (itGroup != t->tree[i].pcount.end()) { pcountSize++; if (pcountSize > 1) { break; } }
223 if (pcountSize == 0) { }
224 else if ((t->tree[i].getBranchLength() != -1) && (pcountSize == 1) && (rootForGrouping[namesOfGroupCombos[h]].count(i) == 0)) { //you have a unique branch length and you are not the root
225 UniqueBL += abs(t->tree[i].getBranchLength());
229 if (pcountSize == 0) { }
230 else if ((t->tree[i].getBranchLength() != -1) && (pcountSize != 0) && (rootForGrouping[namesOfGroupCombos[h]].count(i) == 0)) { //you have a branch length and you are not the root
231 totalBL += abs(t->tree[i].getBranchLength());
234 //cout << UniqueBL << '\t' << totalBL << endl;
235 UW = (UniqueBL / totalBL);
237 if (isnan(UW) || isinf(UW)) { UW = 0; }
244 //if((count % twentyPercent) == 0) { float tempOut = (count / (float)total); if (isnan(tempOut) || isinf(tempOut)) { tempOut = 0.0; } m->mothurOut("Percentage complete: " + toString((int(tempOut) * 100.0))); m->mothurOutEndLine(); }
248 //if((count % twentyPercent) != 0) { float tempOut = (count / (float)total); if (isnan(tempOut) || isinf(tempOut)) { tempOut = 0.0; } m->mothurOut("Percentage complete: " + toString((int(tempOut) * 100.0))); m->mothurOutEndLine(); }
252 catch(exception& e) {
253 m->errorOut(e, "Unweighted", "driver");
257 /**************************************************************************************************/
259 EstOutput Unweighted::getValues(Tree* t, string groupA, string groupB, int p, string o) {
264 TreeMap* tmap = t->getTreeMap();
266 //if the users enters no groups then give them the score of all groups
267 int numGroups = m->getNumGroups();
269 //calculate number of comparsions
271 vector< vector<string> > namesOfGroupCombos;
272 for (int r=0; r<numGroups; r++) {
273 for (int l = 0; l < r; l++) {
275 vector<string> groups; groups.push_back((m->getGroups())[r]); groups.push_back((m->getGroups())[l]);
276 namesOfGroupCombos.push_back(groups);
281 vector<string> groups;
282 if (numGroups == 0) {
283 //get score for all users groups
284 for (int i = 0; i < (tmap->getNamesOfGroups()).size(); i++) {
285 if ((tmap->getNamesOfGroups())[i] != "xxx") {
286 groups.push_back((tmap->getNamesOfGroups())[i]);
289 namesOfGroupCombos.push_back(groups);
291 for (int i = 0; i < m->getNumGroups(); i++) {
292 groups.push_back((m->getGroups())[i]);
294 namesOfGroupCombos.push_back(groups);
298 #if defined (__APPLE__) || (__MACH__) || (linux) || (__linux) || (__linux__) || (__unix__) || (__unix)
300 data = driver(t, namesOfGroupCombos, 0, namesOfGroupCombos.size(), true, tmap);
302 int numPairs = namesOfGroupCombos.size();
304 int numPairsPerProcessor = numPairs / processors;
306 for (int i = 0; i < processors; i++) {
307 int startPos = i * numPairsPerProcessor;
308 if(i == processors - 1){
309 numPairsPerProcessor = numPairs - i * numPairsPerProcessor;
311 lines.push_back(linePair(startPos, numPairsPerProcessor));
314 data = createProcesses(t, namesOfGroupCombos, true, tmap);
319 data = driver(t, namesOfGroupCombos, 0, namesOfGroupCombos.size(), true, tmap);
324 catch(exception& e) {
325 m->errorOut(e, "Unweighted", "getValues");
329 /**************************************************************************************************/
331 EstOutput Unweighted::createProcesses(Tree* t, vector< vector<string> > namesOfGroupCombos, bool usingGroups, TreeMap* tmap) {
333 #if defined (__APPLE__) || (__MACH__) || (linux) || (__linux) || (__linux__) || (__unix__) || (__unix)
335 vector<int> processIDS;
339 //loop through and create all the processes you want
340 while (process != processors) {
344 processIDS.push_back(pid); //create map from line number to pid so you can append files in correct order later
348 myresults = driver(t, namesOfGroupCombos, lines[process].start, lines[process].num, usingGroups, tmap);
350 if (m->control_pressed) { exit(0); }
352 //pass numSeqs to parent
354 string tempFile = outputDir + toString(getpid()) + ".unweighted.results.temp";
355 m->openOutputFile(tempFile, out);
356 out << myresults.size() << endl;
357 for (int i = 0; i < myresults.size(); i++) { out << myresults[i] << '\t'; } out << endl;
362 m->mothurOut("[ERROR]: unable to spawn the necessary processes."); m->mothurOutEndLine();
363 for (int i = 0; i < processIDS.size(); i++) { kill (processIDS[i], SIGINT); }
368 results = driver(t, namesOfGroupCombos, lines[0].start, lines[0].num, usingGroups, tmap);
370 //force parent to wait until all the processes are done
371 for (int i=0;i<(processors-1);i++) {
372 int temp = processIDS[i];
376 if (m->control_pressed) { return results; }
378 //get data created by processes
379 for (int i=0;i<(processors-1);i++) {
381 string s = outputDir + toString(processIDS[i]) + ".unweighted.results.temp";
382 m->openInputFile(s, in);
387 in >> num; m->gobble(in);
389 if (m->control_pressed) { break; }
392 for (int j = 0; j < num; j++) {
395 results.push_back(w);
406 catch(exception& e) {
407 m->errorOut(e, "Unweighted", "createProcesses");
411 /**************************************************************************************************/
412 EstOutput Unweighted::driver(Tree* t, vector< vector<string> > namesOfGroupCombos, int start, int num, bool usingGroups, TreeMap* tmap) {
415 EstOutput results; results.resize(num);
419 Tree* copyTree = new Tree(tmap);
421 for (int h = start; h < (start+num); h++) {
423 if (m->control_pressed) { return results; }
425 //copy random tree passed in
426 copyTree->getCopy(t);
428 //swap labels in the groups you want to compare
429 copyTree->assembleRandomUnifracTree(namesOfGroupCombos[h]);
431 double UniqueBL=0.0000; //a branch length is unique if it's chidren are from the same group
432 double totalBL = 0.00; //all branch lengths
433 double UW = 0.00; //Unweighted Value = UniqueBL / totalBL;
434 //find a node that belongs to one of the groups in this combo
435 int nodeBelonging = -1;
436 for (int g = 0; g < namesOfGroupCombos[h].size(); g++) {
437 if (copyTree->groupNodeInfo[namesOfGroupCombos[h][g]].size() != 0) { nodeBelonging = copyTree->groupNodeInfo[namesOfGroupCombos[h][g]][0]; break; }
441 if (nodeBelonging == -1) {
442 m->mothurOut("[WARNING]: cannot find a nodes in the tree from grouping ");
443 for (int g = 0; g < namesOfGroupCombos[h].size()-1; g++) { m->mothurOut(namesOfGroupCombos[h][g] + "-"); }
444 m->mothurOut(namesOfGroupCombos[h][namesOfGroupCombos[h].size()-1]);
445 m->mothurOut(", skipping."); m->mothurOutEndLine(); results[count] = UW;
448 //if including the root this clears rootForGrouping[namesOfGroupCombos[h]]
449 getRoot(copyTree, nodeBelonging, namesOfGroupCombos[h]);
451 for(int i=0;i<copyTree->getNumNodes();i++){
453 if (m->control_pressed) { return data; }
455 //pcountSize = 0, they are from a branch that is entirely from a group the user doesn't want
456 //pcountSize = 2, not unique to one group
457 //pcountSize = 1, unique to one group
460 for (int j = 0; j < namesOfGroupCombos[h].size(); j++) {
461 map<string, int>::iterator itGroup = copyTree->tree[i].pcount.find(namesOfGroupCombos[h][j]);
462 if (itGroup != copyTree->tree[i].pcount.end()) { pcountSize++; if (pcountSize > 1) { break; } }
466 if (pcountSize == 0) { }
467 else if ((copyTree->tree[i].getBranchLength() != -1) && (pcountSize == 1) && (rootForGrouping[namesOfGroupCombos[h]].count(i) == 0)) { //you have a unique branch length and you are not the root
468 UniqueBL += abs(copyTree->tree[i].getBranchLength());
472 if (pcountSize == 0) { }
473 else if ((copyTree->tree[i].getBranchLength() != -1) && (pcountSize != 0) && (rootForGrouping[namesOfGroupCombos[h]].count(i) == 0)) { //you have a branch length and you are not the root
474 totalBL += abs(copyTree->tree[i].getBranchLength());
478 //cout << UniqueBL << '\t' << totalBL << endl;
479 UW = (UniqueBL / totalBL);
481 if (isnan(UW) || isinf(UW)) { UW = 0; }
493 catch(exception& e) {
494 m->errorOut(e, "Unweighted", "driver");
498 /**************************************************************************************************/
499 int Unweighted::getRoot(Tree* t, int v, vector<string> grouping) {
501 //you are a leaf so get your parent
502 int index = t->tree[v].getParent();
505 rootForGrouping[grouping].clear();
508 //my parent is a potential root
509 rootForGrouping[grouping].insert(index);
511 //while you aren't at root
512 while(t->tree[index].getParent() != -1){
513 //cout << index << endl;
514 if (m->control_pressed) { return 0; }
516 //am I the root for this grouping? if so I want to stop "early"
517 //does my sibling have descendants from the users groups?
518 //if so I am not the root
519 int parent = t->tree[index].getParent();
520 int lc = t->tree[parent].getLChild();
521 int rc = t->tree[parent].getRChild();
524 if (lc == index) { sib = rc; }
526 map<string, int>::iterator itGroup;
528 for (int j = 0; j < grouping.size(); j++) {
529 map<string, int>::iterator itGroup = t->tree[sib].pcount.find(grouping[j]);
530 if (itGroup != t->tree[sib].pcount.end()) { pcountSize++; if (pcountSize > 1) { break; } }
533 //if yes, I am not the root
534 if (pcountSize != 0) {
535 rootForGrouping[grouping].clear();
536 rootForGrouping[grouping].insert(parent);
542 //get all nodes above the root to add so we don't add their u values above
543 index = *(rootForGrouping[grouping].begin());
544 while(t->tree[index].getParent() != -1){
545 int parent = t->tree[index].getParent();
546 rootForGrouping[grouping].insert(parent);
547 //cout << parent << " in root" << endl;
554 catch(exception& e) {
555 m->errorOut(e, "Unweighted", "getRoot");
559 /**************************************************************************************************/