From f099fdc1e3a0d7b75d780a164e5bdb93496a7a1d Mon Sep 17 00:00:00 2001 From: westcott Date: Tue, 31 Aug 2010 16:03:14 +0000 Subject: [PATCH] working on parallelizing unifrac.weighted. --- parsimony.h | 2 +- sharedthetayc.cpp | 4 +- treecalculator.h | 8 +- unifracweightedcommand.cpp | 247 +++++++++++++++++------- unifracweightedcommand.h | 10 + unweighted.cpp | 4 +- weighted.cpp | 371 ++++++++++++++++++++++--------------- weighted.h | 19 +- 8 files changed, 440 insertions(+), 225 deletions(-) diff --git a/parsimony.h b/parsimony.h index d4f71ab..74ebefd 100644 --- a/parsimony.h +++ b/parsimony.h @@ -23,7 +23,7 @@ class Parsimony : public TreeCalculator { Parsimony(TreeMap* t) : tmap(t) {}; ~Parsimony() {}; EstOutput getValues(Tree*); - EstOutput getValues(Tree*, string, string) { return data; }; + //EstOutput getValues(Tree*, string, string) { return data; } private: GlobalData* globaldata; diff --git a/sharedthetayc.cpp b/sharedthetayc.cpp index 4d8a5d9..ed18d92 100644 --- a/sharedthetayc.cpp +++ b/sharedthetayc.cpp @@ -31,8 +31,8 @@ EstOutput ThetaYC::getValues(vector shared) { //get the total values we need to calculate the theta denominator sums for (int i = 0; i < shared[0]->size(); i++) { //store in temps to avoid multiple repetitive function calls - Atotal += (float)shared[0]->getAbundance(i); - Btotal += (float)shared[1]->getAbundance(i); + Atotal += (double)shared[0]->getAbundance(i); + Btotal += (double)shared[1]->getAbundance(i); } //calculate the theta denominator sums diff --git a/treecalculator.h b/treecalculator.h index 86cb527..b2058cb 100644 --- a/treecalculator.h +++ b/treecalculator.h @@ -25,11 +25,15 @@ class TreeCalculator { public: TreeCalculator(){ m = MothurOut::getInstance(); } TreeCalculator(string n) : name(n) {}; + virtual ~TreeCalculator(){}; - virtual EstOutput getValues(Tree*) = 0; - virtual EstOutput getValues(Tree*, string, string) = 0; + virtual EstOutput getValues(Tree*) { return data; } + virtual EstOutput getValues(Tree*, int, string) { return data; } + virtual EstOutput getValues(Tree*, string, string) { return data; } + virtual EstOutput getValues(Tree*, string, string, vector&) { return data; } virtual string getName() { return name; } + virtual vector getBranchLengthSums(Tree*) { vector temp; return temp; } protected: EstOutput data; diff --git a/unifracweightedcommand.cpp b/unifracweightedcommand.cpp index 0201776..65b0a50 100644 --- a/unifracweightedcommand.cpp +++ b/unifracweightedcommand.cpp @@ -21,7 +21,7 @@ UnifracWeightedCommand::UnifracWeightedCommand(string option) { else { //valid paramters for this command - string Array[] = {"groups","iters","distance","random","outputdir","inputdir"}; + string Array[] = {"groups","iters","distance","random","processors","outputdir","inputdir"}; vector myArray (Array, Array+(sizeof(Array)/sizeof(string))); OptionParser parser(option); @@ -61,6 +61,9 @@ UnifracWeightedCommand::UnifracWeightedCommand(string option) { temp = validParameter.validFile(parameters, "random", false); if (temp == "not found") { temp = "F"; } random = m->isTrue(temp); + temp = validParameter.validFile(parameters, "processors", false); if (temp == "not found"){ temp = "1"; } + convert(temp, processors); + if (!random) { iters = 0; } //turn off random calcs @@ -118,9 +121,6 @@ int UnifracWeightedCommand::execute() { int start = time(NULL); - Progress* reading; - if (random) { reading = new Progress("Comparing to random:", iters); } - //get weighted for users tree userData.resize(numComp,0); //data[0] = weightedscore AB, data[1] = weightedscore AC... randomData.resize(numComp,0); //data[0] = weightedscore AB, data[1] = weightedscore AC... @@ -131,13 +131,7 @@ int UnifracWeightedCommand::execute() { //get weighted scores for users trees for (int i = 0; i < T.size(); i++) { - if (m->control_pressed) { - delete randT; - if (random) { delete reading; } - outSum.close(); - for (int i = 0; i < outputNames.size(); i++) { remove(outputNames[i].c_str()); } - return 0; - } + if (m->control_pressed) { delete randT; outSum.close(); for (int i = 0; i < outputNames.size(); i++) { remove(outputNames[i].c_str()); } return 0; } counter = 0; rScores.resize(numComp); //data[0] = weightedscore AB, data[1] = weightedscore AC... @@ -148,11 +142,11 @@ int UnifracWeightedCommand::execute() { outputNames.push_back(outputDir + m->getSimpleName(globaldata->getTreeFile()) + toString(i+1) + ".weighted"); } - userData = weighted->getValues(T[i]); //userData[0] = weightedscore + userData = weighted->getValues(T[i], processors, outputDir); //userData[0] = weightedscore if (m->control_pressed) { delete randT; - if (random) { delete reading; delete output; } + if (random) { delete output; } outSum.close(); for (int i = 0; i < outputNames.size(); i++) { remove(outputNames[i].c_str()); } return 0; @@ -168,51 +162,54 @@ int UnifracWeightedCommand::execute() { utreeScores.push_back(userData[s]); } - //get scores for random trees - for (int j = 0; j < iters; j++) { - int count = 0; - for (int r=0; rgetCopy(T[i]); - - //create a random tree with same topology as T[i], but different labels - randT->assembleRandomUnifracTree(globaldata->Groups[r], globaldata->Groups[l]); - - if (m->control_pressed) { - delete randT; - if (random) { delete reading; delete output; } - outSum.close(); - for (int i = 0; i < outputNames.size(); i++) { remove(outputNames[i].c_str()); } - return 0; + if (random) { + vector sums = weighted->getBranchLengthSums(T[i]); + + //calculate number of comparisons i.e. with groups A,B,C = AB, AC, BC = 3; + vector< vector > namesOfGroupCombos; + for (int a=0; a groups; groups.push_back(globaldata->Groups[a]); groups.push_back(globaldata->Groups[l]); + namesOfGroupCombos.push_back(groups); + } + } + + #if defined (__APPLE__) || (__MACH__) || (linux) || (__linux) + if(processors != 1){ + int numPairs = namesOfGroupCombos.size(); + int numPairsPerProcessor = numPairs / processors; + + for (int i = 0; i < processors; i++) { + int startPos = i * numPairsPerProcessor; + if(i == processors - 1){ + numPairsPerProcessor = numPairs - i * numPairsPerProcessor; + } + lines.push_back(new linePair(startPos, numPairsPerProcessor)); } + } + #endif - - //get wscore of random tree - randomData = weighted->getValues(randT, globaldata->Groups[r], globaldata->Groups[l]); - - if (m->control_pressed) { - delete randT; - if (random) { delete reading; delete output; } - outSum.close(); - for (int i = 0; i < outputNames.size(); i++) { remove(outputNames[i].c_str()); } - return 0; + + //get scores for random trees + for (int j = 0; j < iters; j++) { + int count = 0; + + #if defined (__APPLE__) || (__MACH__) || (linux) || (__linux) + if(processors == 1){ + driver(T[i], randT, namesOfGroupCombos, 0, namesOfGroupCombos.size(), sums, rScores); + }else{ + createProcesses(T[i], randT, namesOfGroupCombos, sums, rScores); + for (int i = 0; i < lines.size(); i++) { delete lines[i]; } lines.clear(); } - - //save scores - rScores[count].push_back(randomData[0]); - count++; - } + #else + driver(T[i], randT, namesOfGroupCombos, 0, namesOfGroupCombos.size(), sums, rScores); + #endif + + if (m->control_pressed) { delete output; outSum.close(); for (int i = 0; i < outputNames.size(); i++) { remove(outputNames[i].c_str()); } return 0; } } - //update progress bar - reading->update(j); - - } - - //removeValidScoresDuplicates(); - //find the signifigance of the score for summary file - if (random) { + + //find the signifigance of the score for summary file for (int f = 0; f < numComp; f++) { //sort random scores sort(rScores[f].begin(), rScores[f].end()); @@ -232,6 +229,7 @@ int UnifracWeightedCommand::execute() { printWeightedFile(); delete output; + } //clear data @@ -241,16 +239,7 @@ int UnifracWeightedCommand::execute() { } - if (m->control_pressed) { - delete randT; - if (random) { delete reading; } - outSum.close(); - for (int i = 0; i < outputNames.size(); i++) { remove(outputNames[i].c_str()); } - return 0; - } - - //finish progress bar - if (random) { reading->finish(); delete reading; } + if (m->control_pressed) { delete randT; outSum.close(); for (int i = 0; i < outputNames.size(); i++) { remove(outputNames[i].c_str()); } return 0; } printWSummaryFile(); @@ -281,6 +270,138 @@ int UnifracWeightedCommand::execute() { exit(1); } } +/**************************************************************************************************/ + +int UnifracWeightedCommand::createProcesses(Tree* t, Tree* randT, vector< vector > namesOfGroupCombos, vector& sums, vector< vector >& scores) { + try { +#if defined (__APPLE__) || (__MACH__) || (linux) || (__linux) + int process = 0; + int num = 0; + vector processIDS; + + EstOutput results; + + //loop through and create all the processes you want + while (process != processors) { + int pid = fork(); + + if (pid > 0) { + processIDS.push_back(pid); //create map from line number to pid so you can append files in correct order later + process++; + }else if (pid == 0){ + driver(t, randT, namesOfGroupCombos, lines[process]->start, lines[process]->num, sums, scores); + + m->mothurOut("Merging results."); m->mothurOutEndLine(); + + //pass numSeqs to parent + ofstream out; + string tempFile = outputDir + toString(getpid()) + ".results.temp"; + m->openOutputFile(tempFile, out); + out << results.size() << endl; + for (int i = lines[process]->start; i < (lines[process]->start + lines[process]->num); i++) { out << results[i] << '\t'; } out << endl; + out.close(); + + exit(0); + }else { m->mothurOut("unable to spawn the necessary processes."); m->mothurOutEndLine(); exit(0); } + } + + //force parent to wait until all the processes are done + for (int i=0;iopenInputFile(s, in); + + vector r; + + //get quantiles + while (!in.eof()) { + int num; + in >> num; + + m->gobble(in); + + double w; + for (int j = 0; j < num; j++) { + in >> w; + r.push_back(w); + } + m->gobble(in); + } + in.close(); + remove(s.c_str()); + + //save quan in quantiles + for (int j = 0; j < r.size(); j++) { + //put all values of r into results + results.push_back(r[j]); + } + } + + m->mothurOut("DONE."); m->mothurOutEndLine(); m->mothurOutEndLine(); + + return 0; +#endif + } + catch(exception& e) { + m->errorOut(e, "UnifracWeightedCommand", "createProcesses"); + exit(1); + } +} + +/**************************************************************************************************/ +int UnifracWeightedCommand::driver(Tree* t, Tree* randT, vector< vector > namesOfGroupCombos, int start, int num, vector& sums, vector< vector >& scores) { + try { + int count = 0; + int total = start+num; + int twentyPercent = (total * 0.20); + + for (int h = start; h < (start+num); h++) { + + if (m->control_pressed) { return 0; } + + //initialize weighted score + string groupA = namesOfGroupCombos[h][0]; + string groupB = namesOfGroupCombos[h][1]; + + //copy T[i]'s info. + randT->getCopy(t); + + //create a random tree with same topology as T[i], but different labels + randT->assembleRandomUnifracTree(groupA, groupB); + + if (m->control_pressed) { delete randT; return 0; } + + + //get wscore of random tree + EstOutput randomData = weighted->getValues(randT, groupA, groupB, sums); + + if (m->control_pressed) { delete randT; return 0; } + + //save scores + scores[h].push_back(randomData[0]); + + count++; + + //report progress + if((h) % twentyPercent == 0){ m->mothurOut("Random comparison percentage complete: " + toString(int((h / (float)total) * 100.0))); m->mothurOutEndLine(); } + } + + m->mothurOut("Random comparison percentage complete: 100"); m->mothurOutEndLine(); + + return 0; + + } + catch(exception& e) { + m->errorOut(e, "UnifracWeightedCommand", "driver"); + exit(1); + } +} /***********************************************************/ void UnifracWeightedCommand::printWeightedFile() { try { diff --git a/unifracweightedcommand.h b/unifracweightedcommand.h index 77dd10d..b6aea7d 100644 --- a/unifracweightedcommand.h +++ b/unifracweightedcommand.h @@ -29,6 +29,13 @@ class UnifracWeightedCommand : public Command { void help(); private: + struct linePair { + int start; + int num; + linePair(int i, int j) : start(i), num(j) {} + }; + vector lines; + GlobalData* globaldata; SharedUtil* util; FileOutput* output; @@ -52,6 +59,7 @@ class UnifracWeightedCommand : public Command { bool abort, phylip, random; string groups, itersString; vector Groups, outputNames; //holds groups to be used + int processors; ofstream outSum; @@ -62,6 +70,8 @@ class UnifracWeightedCommand : public Command { //void removeValidScoresDuplicates(); int findIndex(float, int); void calculateFreqsCumuls(); + int createProcesses(Tree*, Tree*, vector< vector >, vector&, vector< vector >&); + int driver(Tree*, Tree*, vector< vector >, int, int, vector&, vector< vector >&); }; diff --git a/unweighted.cpp b/unweighted.cpp index d9b3b58..c4548c9 100644 --- a/unweighted.cpp +++ b/unweighted.cpp @@ -55,7 +55,7 @@ EstOutput Unweighted::getValues(Tree* t) { int pcountSize = 0; for (int j = 0; j < groups.size(); j++) { map::iterator itGroup = t->tree[i].pcount.find(groups[j]); - if (itGroup != t->tree[i].pcount.end()) { pcountSize++; if (pcountSize > 1) { break; } } + if (itGroup != t->tree[i].pcount.end()) { pcountSize++; } } if (pcountSize == 0) { } @@ -188,7 +188,7 @@ EstOutput Unweighted::getValues(Tree* t, string groupA, string groupB) { int pcountSize = 0; for (int j = 0; j < groups.size(); j++) { map::iterator itGroup = copyTree->tree[i].pcount.find(groups[j]); - if (itGroup != copyTree->tree[i].pcount.end()) { pcountSize++; if (pcountSize > 1) { break; } } + if (itGroup != copyTree->tree[i].pcount.end()) { pcountSize++; } } if (pcountSize == 0) { } diff --git a/weighted.cpp b/weighted.cpp index c8b64e2..d068369 100644 --- a/weighted.cpp +++ b/weighted.cpp @@ -11,11 +11,15 @@ /**************************************************************************************************/ -EstOutput Weighted::getValues(Tree* t) { +EstOutput Weighted::getValues(Tree* t, int p, string o) { try { globaldata = GlobalData::getInstance(); + + data.clear(); //clear out old values int numGroups; vector D; + processors = p; + outputDir = o; numGroups = globaldata->Groups.size(); @@ -24,135 +28,245 @@ EstOutput Weighted::getValues(Tree* t) { if (m->control_pressed) { return data; } //calculate number of comparisons i.e. with groups A,B,C = AB, AC, BC = 3; - int count = 0; + vector< vector > namesOfGroupCombos; for (int i=0; iGroups[i]+globaldata->Groups[l]] = 0.0; - + //WScore[globaldata->Groups[i]+globaldata->Groups[l]] = 0.0; vector groups; groups.push_back(globaldata->Groups[i]); groups.push_back(globaldata->Groups[l]); + namesOfGroupCombos.push_back(groups); + } + } + + #if defined (__APPLE__) || (__MACH__) || (linux) || (__linux) + if(processors == 1){ + data = driver(t, namesOfGroupCombos, 0, namesOfGroupCombos.size(), sums); + }else{ + int numPairs = namesOfGroupCombos.size(); - D.push_back(0.0000); //initialize a spot in D for each combination - - //adding the wieghted sums from group i - for (int j = 0; j < t->groupNodeInfo[groups[0]].size(); j++) { //the leaf nodes that have seqs from group i - map::iterator it = t->tree[t->groupNodeInfo[groups[0]][j]].pcount.find(groups[0]); - int numSeqsInGroupI = it->second; - - double weightedSum = ((numSeqsInGroupI * sums[t->groupNodeInfo[groups[0]][j]]) / (double)tmap->seqsPerGroup[groups[0]]); + int numPairsPerProcessor = numPairs / processors; - D[count] += weightedSum; + for (int i = 0; i < processors; i++) { + int startPos = i * numPairsPerProcessor; + if(i == processors - 1){ + numPairsPerProcessor = numPairs - i * numPairsPerProcessor; + } + lines.push_back(new linePair(startPos, numPairsPerProcessor)); } + + data = createProcesses(t, namesOfGroupCombos, sums); - //adding the wieghted sums from group l - for (int j = 0; j < t->groupNodeInfo[groups[1]].size(); j++) { //the leaf nodes that have seqs from group l - map::iterator it = t->tree[t->groupNodeInfo[groups[1]][j]].pcount.find(groups[1]); - int numSeqsInGroupL = it->second; - - double weightedSum = ((numSeqsInGroupL * sums[t->groupNodeInfo[groups[1]][j]]) / (double)tmap->seqsPerGroup[groups[1]]); - - D[count] += weightedSum; - } + for (int i = 0; i < lines.size(); i++) { delete lines[i]; } lines.clear(); + } + #else + data = driver(t, namesOfGroupCombos, 0, namesOfGroupCombos.size(), sums); + #endif + + return data; + } + catch(exception& e) { + m->errorOut(e, "Weighted", "getValues"); + exit(1); + } +} +/**************************************************************************************************/ - /******************************************************** - //calculate a D value for each group combo - for(int v=0;vgetNumLeaves();v++){ +EstOutput Weighted::createProcesses(Tree* t, vector< vector > namesOfGroupCombos, vector& sums) { + try { +#if defined (__APPLE__) || (__MACH__) || (linux) || (__linux) + int process = 0; + int num = 0; + vector processIDS; + + EstOutput results; + + //loop through and create all the processes you want + while (process != processors) { + int pid = fork(); + + if (pid > 0) { + processIDS.push_back(pid); //create map from line number to pid so you can append files in correct order later + process++; + }else if (pid == 0){ + results = driver(t, namesOfGroupCombos, lines[process]->start, lines[process]->num, sums); - if (m->control_pressed) { return data; } - - //is this sum from a sequence which is in one of the users groups - if (m->inUsersGroups(t->tree[v].getGroup(), globaldata->Groups) == true) { - //is this sum from a sequence which is in this groupCombo - if (m->inUsersGroups(t->tree[v].getGroup(), groups)) { - int numSeqsInGroupI, numSeqsInGroupL; - - map::iterator it; - it = t->tree[v].pcount.find(groups[0]); - if (it != t->tree[v].pcount.end()) { //this leaf node contains seqs from group i - numSeqsInGroupI = it->second; - }else{ numSeqsInGroupI = 0; } - - it = t->tree[v].pcount.find(groups[1]); - if (it != t->tree[v].pcount.end()) { //this leaf node contains seqs from group l - numSeqsInGroupL = it->second; - }else{ numSeqsInGroupL = 0; } - - double weightedSum = ((numSeqsInGroupI * sum) / (double)tmap->seqsPerGroup[groups[0]]) + ((numSeqsInGroupL * sum) / (double)tmap->seqsPerGroup[groups[1]]); + if (m->control_pressed) { exit(0); } + + m->mothurOut("Merging results."); m->mothurOutEndLine(); + + //pass numSeqs to parent + ofstream out; + string tempFile = outputDir + toString(getpid()) + ".results.temp"; + m->openOutputFile(tempFile, out); + out << results.size() << endl; + for (int i = 0; i < results.size(); i++) { out << results[i] << '\t'; } out << endl; + out.close(); + + exit(0); + }else { m->mothurOut("unable to spawn the necessary processes."); m->mothurOutEndLine(); exit(0); } + } + + //force parent to wait until all the processes are done + for (int i=0;icontrol_pressed) { return results; } + + //get data created by processes + for (int i=0;iopenInputFile(s, in); + + vector r; + + //get quantiles + while (!in.eof()) { + int num; + in >> num; + + if (m->control_pressed) { break; } + + m->gobble(in); - //sum /= (double)tmap->seqsPerGroup[t->tree[v].getGroup()]; - - D[count] += weightedSum; - } - } + double w; + for (int j = 0; j < num; j++) { + in >> w; + r.push_back(w); } - /*********************************************************/ - count++; + m->gobble(in); + } + in.close(); + remove(s.c_str()); + + //save quan in quantiles + for (int j = 0; j < r.size(); j++) { + //put all values of r into results + results.push_back(r[j]); } } - data.clear(); //clear out old values - + m->mothurOut("DONE."); m->mothurOutEndLine(); m->mothurOutEndLine(); + + return results; +#endif + } + catch(exception& e) { + m->errorOut(e, "Weighted", "createProcesses"); + exit(1); + } +} +/**************************************************************************************************/ +EstOutput Weighted::driver(Tree* t, vector< vector > namesOfGroupCombos, int start, int num, vector& sums) { + try { + globaldata = GlobalData::getInstance(); + + EstOutput results; + vector D; + + int count = 0; + int total = start+num; + int twentyPercent = (total * 0.20); + + for (int h = start; h < (start+num); h++) { + + if (m->control_pressed) { return results; } + + //initialize weighted score + string groupA = namesOfGroupCombos[h][0]; + string groupB = namesOfGroupCombos[h][1]; + + WScore[groupA+groupB] = 0.0; + D.push_back(0.0000); //initialize a spot in D for each combination + + //adding the wieghted sums from group i + for (int j = 0; j < t->groupNodeInfo[groupA].size(); j++) { //the leaf nodes that have seqs from group i + map::iterator it = t->tree[t->groupNodeInfo[groupA][j]].pcount.find(groupA); + int numSeqsInGroupI = it->second; + + double weightedSum = ((numSeqsInGroupI * sums[t->groupNodeInfo[groupA][j]]) / (double)tmap->seqsPerGroup[groupA]); + + D[count] += weightedSum; + } + + //adding the wieghted sums from group l + for (int j = 0; j < t->groupNodeInfo[groupB].size(); j++) { //the leaf nodes that have seqs from group l + map::iterator it = t->tree[t->groupNodeInfo[groupB][j]].pcount.find(groupB); + int numSeqsInGroupL = it->second; + + double weightedSum = ((numSeqsInGroupL * sums[t->groupNodeInfo[groupB][j]]) / (double)tmap->seqsPerGroup[groupB]); + + D[count] += weightedSum; + } + count++; + + //report progress + if((h) % twentyPercent == 0){ m->mothurOut("Percentage complete: " + toString(int((h / (float)total) * 100.0))); m->mothurOutEndLine(); } + } + + m->mothurOut("Percentage complete: 100"); m->mothurOutEndLine(); + + //calculate u for the group comb for(int i=0;igetNumNodes();i++){ - //calculate weighted score for each of the group comb i.e. with groups A,B,C = AB, AC, BC. - for (int b=0; bcontrol_pressed) { return data; } - - //calculate a u value for each combo - double u; - //does this node have descendants from group b-1 - it = t->tree[i].pcount.find(globaldata->Groups[b]); - //if it does u = # of its descendants with a certain group / total number in tree with a certain group - if (it != t->tree[i].pcount.end()) { - u = (double) t->tree[i].pcount[globaldata->Groups[b]] / (double) tmap->seqsPerGroup[globaldata->Groups[b]]; - }else { u = 0.00; } - - //does this node have descendants from group l - it = t->tree[i].pcount.find(globaldata->Groups[l]); - //if it does subtract their percentage from u - if (it != t->tree[i].pcount.end()) { - u -= (double) t->tree[i].pcount[globaldata->Groups[l]] / (double) tmap->seqsPerGroup[globaldata->Groups[l]]; - } - - u = abs(u * t->tree[i].getBranchLength()); - - //save groupcombs u value - WScore[globaldata->Groups[b]+globaldata->Groups[l]] += u; - /*********************************************************/ + string groupA = namesOfGroupCombos[h][0]; + string groupB = namesOfGroupCombos[h][1]; + + if (m->control_pressed) { return results; } + + double u; + //does this node have descendants from groupA + it = t->tree[i].pcount.find(groupA); + //if it does u = # of its descendants with a certain group / total number in tree with a certain group + if (it != t->tree[i].pcount.end()) { + u = (double) t->tree[i].pcount[groupA] / (double) tmap->seqsPerGroup[groupA]; + }else { u = 0.00; } + + + //does this node have descendants from group l + it = t->tree[i].pcount.find(groupB); + //if it does subtract their percentage from u + if (it != t->tree[i].pcount.end()) { + u -= (double) t->tree[i].pcount[groupB] / (double) tmap->seqsPerGroup[groupB]; } + + u = abs(u * t->tree[i].getBranchLength()); + + //save groupcombs u value + WScore[(groupA+groupB)] += u; } } - - //calculate weighted score for each group combination + + /********************************************************/ + + //calculate weighted score for the group combination double UN; count = 0; - for (int i=0; iGroups[i]+globaldata->Groups[l]] / D[count]); + for (int h = start; h < (start+num); h++) { + UN = (WScore[namesOfGroupCombos[h][0]+namesOfGroupCombos[h][1]] / D[count]); - if (isnan(UN) || isinf(UN)) { UN = 0; } - data.push_back(UN); - count++; - } + if (isnan(UN) || isinf(UN)) { UN = 0; } + results.push_back(UN); + count++; } - return data; + + return results; } catch(exception& e) { - m->errorOut(e, "Weighted", "getValues"); + m->errorOut(e, "Weighted", "driver"); exit(1); } } - /**************************************************************************************************/ -EstOutput Weighted::getValues(Tree* t, string groupA, string groupB) { +EstOutput Weighted::getValues(Tree* t, string groupA, string groupB, vector& sums) { try { globaldata = GlobalData::getInstance(); data.clear(); //clear out old values - vector sums = getBranchLengthSums(t); - if (m->control_pressed) { return data; } //initialize weighted score @@ -161,72 +275,25 @@ EstOutput Weighted::getValues(Tree* t, string groupA, string groupB) { vector groups; groups.push_back(groupA); groups.push_back(groupB); - //adding the wieghted sums from groupA - for (int j = 0; j < t->groupNodeInfo[groupA].size(); j++) { //the leaf nodes that have seqs from group i - map::iterator it = t->tree[j].pcount.find(groupA); - int numSeqsInGroupA = it->second; + //adding the wieghted sums from group i + for (int j = 0; j < t->groupNodeInfo[groups[0]].size(); j++) { //the leaf nodes that have seqs from group i + map::iterator it = t->tree[t->groupNodeInfo[groups[0]][j]].pcount.find(groups[0]); + int numSeqsInGroupI = it->second; - double weightedSum = ((numSeqsInGroupA * sums[j]) / (double)tmap->seqsPerGroup[groupA]); + double weightedSum = ((numSeqsInGroupI * sums[t->groupNodeInfo[groups[0]][j]]) / (double)tmap->seqsPerGroup[groups[0]]); D += weightedSum; } - //adding the wieghted sums from groupB - for (int j = 0; j < t->groupNodeInfo[groupB].size(); j++) { //the leaf nodes that have seqs from group l - map::iterator it = t->tree[j].pcount.find(groupB); - int numSeqsInGroupB = it->second; + //adding the wieghted sums from group l + for (int j = 0; j < t->groupNodeInfo[groups[1]].size(); j++) { //the leaf nodes that have seqs from group l + map::iterator it = t->tree[t->groupNodeInfo[groups[1]][j]].pcount.find(groups[1]); + int numSeqsInGroupL = it->second; - double weightedSum = ((numSeqsInGroupB * sums[j]) / (double)tmap->seqsPerGroup[groupB]); + double weightedSum = ((numSeqsInGroupL * sums[t->groupNodeInfo[groups[1]][j]]) / (double)tmap->seqsPerGroup[groups[1]]); D += weightedSum; } - - - /******************************************************** - //calculate a D value for the group combo - for(int v=0;vgetNumLeaves();v++){ - if (m->control_pressed) { return data; } - - int index = v; - double sum = 0.0000; - - //while you aren't at root - while(t->tree[index].getParent() != -1){ - - //if you have a BL - if(t->tree[index].getBranchLength() != -1){ - sum += abs(t->tree[index].getBranchLength()); - } - index = t->tree[index].getParent(); - } - - //get last breanch length added - if(t->tree[index].getBranchLength() != -1){ - sum += abs(t->tree[index].getBranchLength()); - } - - if (m->inUsersGroups(t->tree[v].getGroup(), groups)) { - int numSeqsInGroupI, numSeqsInGroupL; - - map::iterator it; - it = t->tree[v].pcount.find(groups[0]); - if (it != t->tree[v].pcount.end()) { //this leaf node contains seqs from group i - numSeqsInGroupI = it->second; - }else{ numSeqsInGroupI = 0; } - - it = t->tree[v].pcount.find(groups[1]); - if (it != t->tree[v].pcount.end()) { //this leaf node contains seqs from group l - numSeqsInGroupL = it->second; - }else{ numSeqsInGroupL = 0; } - - double weightedSum = ((numSeqsInGroupI * sum) / (double)tmap->seqsPerGroup[groups[0]]) + ((numSeqsInGroupL * sum) / (double)tmap->seqsPerGroup[groups[1]]); - - //sum /= (double)tmap->seqsPerGroup[t->tree[v].getGroup()]; - - D += weightedSum; - } - } - /********************************************************/ //calculate u for the group comb for(int i=0;igetNumNodes();i++){ diff --git a/weighted.h b/weighted.h index fedfa0c..8a8ad8f 100644 --- a/weighted.h +++ b/weighted.h @@ -21,17 +21,30 @@ class Weighted : public TreeCalculator { public: Weighted(TreeMap* t) : tmap(t) {}; ~Weighted() {}; - EstOutput getValues(Tree*); - EstOutput getValues(Tree*, string, string); + + EstOutput getValues(Tree*, string, string, vector&); + EstOutput getValues(Tree*, int, string); + + vector getBranchLengthSums(Tree*); private: + struct linePair { + int start; + int num; + linePair(int i, int j) : start(i), num(j) {} + }; + vector lines; + GlobalData* globaldata; EstOutput data; TreeMap* tmap; map::iterator it; map WScore; //a score for each group combination i.e. AB, AC, BC. + int processors; + string outputDir; - vector getBranchLengthSums(Tree*); + EstOutput driver(Tree*, vector< vector >, int, int, vector&); + EstOutput createProcesses(Tree*, vector< vector >, vector&); }; /***********************************************************************/ -- 2.39.2