indicatorValues = getValues(groupings, indicatorGroups, randomGroupingsMap);
- pValues = getPValues(groupings, randomGroupingsMap, lookup.size(), indicatorValues);
+ pValues = getPValues(groupings, lookup.size(), indicatorValues);
}else {
vector< vector<SharedRAbundFloatVector*> > groupings;
set<string> groupsAlreadyAdded;
indicatorValues = getValues(groupings, indicatorGroups, randomGroupingsMap);
- pValues = getPValues(groupings, randomGroupingsMap, lookupFloat.size(), indicatorValues);
+ pValues = getPValues(groupings, lookupFloat.size(), indicatorValues);
}
if (m->control_pressed) { out.close(); return 0; }
indicatorValues = getValues(groupings, indicatorGroups, randomGroupingsMap);
- pValues = getPValues(groupings, randomGroupingsMap, lookup.size(), indicatorValues);
+ pValues = getPValues(groupings, lookup.size(), indicatorValues);
}else {
vector< vector<SharedRAbundFloatVector*> > groupings;
indicatorValues = getValues(groupings, indicatorGroups, randomGroupingsMap);
- pValues = getPValues(groupings, randomGroupingsMap, lookupFloat.size(), indicatorValues);
+ pValues = getPValues(groupings, lookupFloat.size(), indicatorValues);
}
if (m->control_pressed) { out.close(); return 0; }
}
}
//**********************************************************************************************************************
-vector<float> IndicatorCommand::driver(vector< vector<SharedRAbundFloatVector*> >& groupings, map< vector<int>, vector<int> > groupingsMap, int num, vector<float> indicatorValues, int numIters){
+vector<float> IndicatorCommand::driver(vector< vector<SharedRAbundFloatVector*> >& groupings, int num, vector<float> indicatorValues, int numIters){
try {
vector<float> pvalues;
pvalues.resize(indicatorValues.size(), 0);
for(int i=0;i<numIters;i++){
if (m->control_pressed) { break; }
- groupingsMap = randomizeGroupings(groupings, num);
+ map< vector<int>, vector<int> > groupingsMap = randomizeGroupings(groupings, num);
vector<float> randomIndicatorValues = getValues(groupings, notUsedGroupings, groupingsMap);
for (int j = 0; j < indicatorValues.size(); j++) {
}
}
//**********************************************************************************************************************
-vector<float> IndicatorCommand::getPValues(vector< vector<SharedRAbundFloatVector*> >& groupings, map< vector<int>, vector<int> > groupingsMap, int num, vector<float> indicatorValues){
+vector<float> IndicatorCommand::getPValues(vector< vector<SharedRAbundFloatVector*> >& groupings, int num, vector<float> indicatorValues){
try {
vector<float> pvalues;
-
-#if defined (__APPLE__) || (__MACH__) || (linux) || (__linux) || (__linux__) || (__unix__) || (__unix)
+
if(processors == 1){
- pvalues = driver(groupings, groupingsMap, num, indicatorValues, iters);
- for (int i = 0; i < pvalues.size(); i++) { pvalues[i] /= (double)iters; }
+ pvalues = driver(groupings, num, indicatorValues, iters);
+ for (int i = 0; i < pvalues.size(); i++) { pvalues[i] /= (double)iters; }
}else{
+ //divide iters between processors
+ vector<int> procIters;
+ int numItersPerProcessor = iters / processors;
+
+ //divide iters between processes
+ for (int h = 0; h < processors; h++) {
+ if(h == processors - 1){ numItersPerProcessor = iters - h * numItersPerProcessor; }
+ procIters.push_back(numItersPerProcessor);
+ }
+
+ vector<int> processIDS;
+ int process = 1;
- //divide iters between processors
- int numItersPerProcessor = iters / processors;
-
- vector<int> processIDS;
- int process = 1;
- int num = 0;
-
+#if defined (__APPLE__) || (__MACH__) || (linux) || (__linux) || (__linux__) || (__unix__) || (__unix)
+
//loop through and create all the processes you want
while (process != processors) {
int pid = fork();
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){
- pvalues = driver(groupings, groupingsMap, num, indicatorValues, numItersPerProcessor);
+ pvalues = driver(groupings, num, indicatorValues, procIters[process]);
//pass pvalues to parent
ofstream out;
}
//do my part
- //special case for last processor in case it doesn't divide evenly
- numItersPerProcessor = iters - ((processors-1) * numItersPerProcessor);
- pvalues = driver(groupings, groupingsMap, num, indicatorValues, numItersPerProcessor);
+ pvalues = driver(groupings, num, indicatorValues, procIters[0]);
//force parent to wait until all the processes are done
for (int i=0;i<processIDS.size();i++) {
}
in.close(); m->mothurRemove(tempFile);
}
- for (int i = 0; i < pvalues.size(); i++) { pvalues[i] /= (double)iters; }
- }
+ for (int i = 0; i < pvalues.size(); i++) { pvalues[i] /= (double)iters; }
#else
- pvalues = driver(groupings, groupingsMap, num, indicatorValues, iters);
- for (int i = 0; i < pvalues.size(); i++) { pvalues[i] /= (double)iters; }
+
+ //fill in functions
+ vector<indicatorData*> pDataArray;
+ DWORD dwThreadIdArray[processors-1];
+ HANDLE hThreadArray[processors-1];
+
+ //Create processor worker threads.
+ for( int i=1; i<processors; i++ ){
+
+ //make copy of lookup so we don't get access violations
+ vector< vector<SharedRAbundFloatVector*> > newGroupings;
+
+ for (int k = 0; k < groupings.size(); k++) {
+ vector<SharedRAbundFloatVector*> newLookup;
+ for (int l = 0; l < groupings[k].size(); l++) {
+ SharedRAbundFloatVector* temp = new SharedRAbundFloatVector();
+ temp->setLabel(groupings[k][l]->getLabel());
+ temp->setGroup(groupings[k][l]->getGroup());
+ newLookup.push_back(temp);
+ }
+ newGroupings.push_back(newLookup);
+ }
+
+ //for each bin
+ for (int l = 0; l < groupings.size(); l++) {
+ for (int k = 0; k < groupings[l][0]->getNumBins(); k++) {
+ if (m->control_pressed) { for (int j = 0; j < newGroupings.size(); j++) { for (int u = 0; u < newGroupings[j].size(); u++) { delete newGroupings[j][u]; } } return pvalues; }
+
+ for (int j = 0; j < groupings[l].size(); j++) { newGroupings[l][j]->push_back(groupings[l][j]->getAbundance(k), groupings[l][j]->getGroup()); }
+ }
+ }
+
+ vector<float> copyIValues = indicatorValues;
+
+ indicatorData* temp = new indicatorData(m, procIters[i], newGroupings, num, copyIValues);
+ pDataArray.push_back(temp);
+ processIDS.push_back(i);
+
+ hThreadArray[i-1] = CreateThread(NULL, 0, MyIndicatorThreadFunction, pDataArray[i-1], 0, &dwThreadIdArray[i-1]);
+ }
+
+ //do my part
+ pvalues = driver(groupings, num, indicatorValues, procIters[0]);
+
+ //Wait until all threads have terminated.
+ WaitForMultipleObjects(processors-1, hThreadArray, TRUE, INFINITE);
+
+ //Close all thread handles and free memory allocations.
+ for(int i=0; i < pDataArray.size(); i++){
+ for (int j = 0; j < pDataArray[i]->pvalues.size(); j++) { pvalues[j] += pDataArray[i]->pvalues[j]; }
+
+ for (int l = 0; l < pDataArray[i]->groupings.size(); l++) {
+ for (int j = 0; j < pDataArray[i]->groupings[l].size(); j++) { delete pDataArray[i]->groupings[l][j]; }
+ }
+
+ CloseHandle(hThreadArray[i]);
+ delete pDataArray[i];
+ }
+
+ for (int i = 0; i < pvalues.size(); i++) { pvalues[i] /= (double)iters; }
#endif
+ }
+
return pvalues;
}
//**********************************************************************************************************************
//same as above, just data type difference
-vector<float> IndicatorCommand::driver(vector< vector<SharedRAbundVector*> >& groupings, map< vector<int>, vector<int> > groupingsMap, int num, vector<float> indicatorValues, int numIters){
+vector<float> IndicatorCommand::driver(vector< vector<SharedRAbundVector*> >& groupings, int num, vector<float> indicatorValues, int numIters){
try {
vector<float> pvalues;
pvalues.resize(indicatorValues.size(), 0);
for(int i=0;i<numIters;i++){
if (m->control_pressed) { break; }
- groupingsMap = randomizeGroupings(groupings, num);
+ map< vector<int>, vector<int> > groupingsMap = randomizeGroupings(groupings, num);
vector<float> randomIndicatorValues = getValues(groupings, notUsedGroupings, groupingsMap);
for (int j = 0; j < indicatorValues.size(); j++) {
}
//**********************************************************************************************************************
//same as above, just data type difference
-vector<float> IndicatorCommand::getPValues(vector< vector<SharedRAbundVector*> >& groupings, map< vector<int>, vector<int> > groupingsMap, int num, vector<float> indicatorValues){
+vector<float> IndicatorCommand::getPValues(vector< vector<SharedRAbundVector*> >& groupings, int num, vector<float> indicatorValues){
try {
vector<float> pvalues;
-
-#if defined (__APPLE__) || (__MACH__) || (linux) || (__linux) || (__linux__) || (__unix__) || (__unix)
+
if(processors == 1){
- pvalues = driver(groupings, groupingsMap, num, indicatorValues, iters);
- for (int i = 0; i < pvalues.size(); i++) { pvalues[i] /= (double)iters; }
+ pvalues = driver(groupings, num, indicatorValues, iters);
+ for (int i = 0; i < pvalues.size(); i++) { pvalues[i] /= (double)iters; }
}else{
+ //divide iters between processors
+ vector<int> procIters;
+ int numItersPerProcessor = iters / processors;
+
+ //divide iters between processes
+ for (int h = 0; h < processors; h++) {
+ if(h == processors - 1){ numItersPerProcessor = iters - h * numItersPerProcessor; }
+ procIters.push_back(numItersPerProcessor);
+ }
+
+ vector<int> processIDS;
+ int process = 1;
- //divide iters between processors
- int numItersPerProcessor = iters / processors;
-
- vector<int> processIDS;
- int process = 1;
-
+#if defined (__APPLE__) || (__MACH__) || (linux) || (__linux) || (__linux__) || (__unix__) || (__unix)
+
//loop through and create all the processes you want
while (process != processors) {
int pid = fork();
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){
- pvalues = driver(groupings, groupingsMap, num, indicatorValues, numItersPerProcessor);
+ pvalues = driver(groupings, num, indicatorValues, procIters[process]);
//pass pvalues to parent
ofstream out;
out.close();
exit(0);
- }else {
- m->mothurOut("[ERROR]: unable to spawn the necessary processes."); m->mothurOutEndLine();
+ }else {
+ m->mothurOut("[ERROR]: unable to spawn the necessary processes."); m->mothurOutEndLine();
for (int i = 0; i < processIDS.size(); i++) { kill (processIDS[i], SIGINT); }
exit(0);
}
}
//do my part
- //special case for last processor in case it doesn't divide evenly
- numItersPerProcessor = iters - ((processors-1) * numItersPerProcessor);
- pvalues = driver(groupings, groupingsMap, num, indicatorValues, numItersPerProcessor);
+ pvalues = driver(groupings, num, indicatorValues, procIters[0]);
//force parent to wait until all the processes are done
- for (int i=0;i<processIDS.size();i++) {
+ for (int i=0;i<processIDS.size();i++) {
int temp = processIDS[i];
wait(&temp);
}
- //combine results
+ //combine results
for (int i = 0; i < processIDS.size(); i++) {
ifstream in;
string tempFile = toString(processIDS[i]) + ".pvalues.temp";
m->openInputFile(tempFile, in);
////// to do ///////////
- int numTemp; numTemp = 0;
+ int numTemp; numTemp = 0;
for (int j = 0; j < pvalues.size(); j++) {
in >> numTemp; m->gobble(in);
pvalues[j] += numTemp;
}
in.close(); m->mothurRemove(tempFile);
}
- for (int i = 0; i < pvalues.size(); i++) { pvalues[i] /= (double)iters; }
- }
+ for (int i = 0; i < pvalues.size(); i++) { pvalues[i] /= (double)iters; }
#else
- pvalues = driver(groupings, groupingsMap, num, indicatorValues, iters);
- for (int i = 0; i < pvalues.size(); i++) { pvalues[i] /= (double)iters; }
+
+ //fill in functions
+ vector<indicatorData*> pDataArray;
+ DWORD dwThreadIdArray[processors-1];
+ HANDLE hThreadArray[processors-1];
+
+ //Create processor worker threads.
+ for( int i=1; i<processors; i++ ){
+
+ //make copy of lookup so we don't get access violations
+ vector< vector<SharedRAbundFloatVector*> > newGroupings;
+
+ for (int k = 0; k < groupings.size(); k++) {
+ vector<SharedRAbundFloatVector*> newLookup;
+ for (int l = 0; l < groupings[k].size(); l++) {
+ SharedRAbundFloatVector* temp = new SharedRAbundFloatVector();
+ temp->setLabel(groupings[k][l]->getLabel());
+ temp->setGroup(groupings[k][l]->getGroup());
+ newLookup.push_back(temp);
+ }
+ newGroupings.push_back(newLookup);
+ }
+
+ //for each bin
+ for (int l = 0; l < groupings.size(); l++) {
+ for (int k = 0; k < groupings[l][0]->getNumBins(); k++) {
+ if (m->control_pressed) { for (int j = 0; j < newGroupings.size(); j++) { for (int u = 0; u < newGroupings[j].size(); u++) { delete newGroupings[j][u]; } } return pvalues; }
+
+ for (int j = 0; j < groupings[l].size(); j++) { newGroupings[l][j]->push_back((float)(groupings[l][j]->getAbundance(k)), groupings[l][j]->getGroup()); }
+ }
+ }
+
+ vector<float> copyIValues = indicatorValues;
+
+ indicatorData* temp = new indicatorData(m, procIters[i], newGroupings, num, copyIValues);
+ pDataArray.push_back(temp);
+ processIDS.push_back(i);
+
+ hThreadArray[i-1] = CreateThread(NULL, 0, MyIndicatorThreadFunction, pDataArray[i-1], 0, &dwThreadIdArray[i-1]);
+ }
+
+ //do my part
+ pvalues = driver(groupings, num, indicatorValues, procIters[0]);
+
+ //Wait until all threads have terminated.
+ WaitForMultipleObjects(processors-1, hThreadArray, TRUE, INFINITE);
+
+ //Close all thread handles and free memory allocations.
+ for(int i=0; i < pDataArray.size(); i++){
+ for (int j = 0; j < pDataArray[i]->pvalues.size(); j++) { pvalues[j] += pDataArray[i]->pvalues[j]; }
+
+ for (int l = 0; l < pDataArray[i]->groupings.size(); l++) {
+ for (int j = 0; j < pDataArray[i]->groupings[l].size(); j++) { delete pDataArray[i]->groupings[l][j]; }
+ }
+
+ CloseHandle(hThreadArray[i]);
+ delete pDataArray[i];
+ }
+
+ for (int i = 0; i < pvalues.size(); i++) { pvalues[i] /= (double)iters; }
#endif
+ }
+
return pvalues;
}
catch(exception& e) {
- m->errorOut(e, "IndicatorCommand", "getPValues");
+ m->errorOut(e, "IndicatorCommand", "getPValues");
exit(1);
}
}
map< vector<int>, vector<int> > randomizeGroupings(vector< vector<SharedRAbundVector*> >&, int);
map< vector<int>, vector<int> > randomizeGroupings(vector< vector<SharedRAbundFloatVector*> >&, int);
- vector<float> driver(vector< vector<SharedRAbundFloatVector*> >&, map< vector<int>, vector<int> >, int, vector<float>, int);
- vector<float> driver(vector< vector<SharedRAbundVector*> >&, map< vector<int>, vector<int> >, int, vector<float>, int);
+ vector<float> driver(vector< vector<SharedRAbundFloatVector*> >&, int, vector<float>, int);
+ vector<float> driver(vector< vector<SharedRAbundVector*> >&, int, vector<float>, int);
- vector<float> getPValues(vector< vector<SharedRAbundFloatVector*> >&, map< vector<int>, vector<int> >, int, vector<float>);
- vector<float> getPValues(vector< vector<SharedRAbundVector*> >&, map< vector<int>, vector<int> >, int, vector<float>);
+ vector<float> getPValues(vector< vector<SharedRAbundFloatVector*> >&, int, vector<float>);
+ vector<float> getPValues(vector< vector<SharedRAbundVector*> >&, int, vector<float>);
};
+/**************************************************************************************************/
+
+struct indicatorData {
+ vector< vector<SharedRAbundFloatVector*> > groupings;
+ MothurOut* m;
+ int iters, num;
+ vector<float> indicatorValues;
+ vector<float> pvalues;
+
+ indicatorData(){}
+ indicatorData(MothurOut* mout, int it, vector< vector<SharedRAbundFloatVector*> > ng, int n, vector<float> iv) {
+ m = mout;
+ iters = it;
+ groupings = ng;
+ indicatorValues = iv;
+ num = n;
+ }
+};
+/**************************************************************************************************/
+#if defined (__APPLE__) || (__MACH__) || (linux) || (__linux) || (__linux__) || (__unix__) || (__unix)
+#else
+static DWORD WINAPI MyIndicatorThreadFunction(LPVOID lpParam){
+ indicatorData* pDataArray;
+ pDataArray = (indicatorData*)lpParam;
+
+ try {
+
+ pDataArray->pvalues.resize(pDataArray->indicatorValues.size(), 0);
+
+ for(int i=0;i<pDataArray->iters;i++){
+ if (pDataArray->m->control_pressed) { break; }
+
+ //groupingsMap = randomizeGroupings(groupings, num);
+ ///////////////////////////////////////////////////////////////////////
+ map< vector<int>, vector<int> > randomGroupings;
+
+ for (int j = 0; j < pDataArray->num; j++) {
+
+ //get random groups to swap to switch with
+ //generate random int between 0 and groupings.size()-1
+ int z = pDataArray->m->getRandomIndex(pDataArray->groupings.size()-1);
+ int x = pDataArray->m->getRandomIndex(pDataArray->groupings.size()-1);
+ int a = pDataArray->m->getRandomIndex(pDataArray->groupings[z].size()-1);
+ int b = pDataArray->m->getRandomIndex(pDataArray->groupings[x].size()-1);
+ //cout << i << '\t' << z << '\t' << x << '\t' << a << '\t' << b << endl;
+
+ vector<int> from;
+ vector<int> to;
+
+ from.push_back(z); from.push_back(a);
+ to.push_back(x); to.push_back(b);
+
+ randomGroupings[from] = to;
+ }
+ ///////////////////////////////////////////////////////////////////////
+
+ //vector<float> randomIndicatorValues = getValues(groupings, notUsedGroupings, randomGroupings);
+ ///////////////////////////////////////////////////////////////////////
+ vector<float> randomIndicatorValues;
+ map< vector<int>, vector<int> >::iterator it;
+
+ //for each otu
+ for (int i = 0; i < pDataArray->groupings[0][0]->getNumBins(); i++) {
+
+ if (pDataArray->m->control_pressed) { return 0; }
+
+ vector<float> terms;
+ float AijDenominator = 0.0;
+ vector<float> Bij;
+
+ //get overall abundance of each grouping
+ for (int j = 0; j < pDataArray->groupings.size(); j++) {
+
+ float totalAbund = 0;
+ int numNotZero = 0;
+ for (int k = 0; k < pDataArray->groupings[j].size(); k++) {
+ vector<int> temp; temp.push_back(j); temp.push_back(k);
+ it = randomGroupings.find(temp);
+
+ if (it == randomGroupings.end()) { //this one didnt get moved
+ totalAbund += pDataArray->groupings[j][k]->getAbundance(i);
+ if (pDataArray->groupings[j][k]->getAbundance(i) != 0.0) { numNotZero++; }
+ }else {
+ totalAbund += pDataArray->groupings[(it->second)[0]][(it->second)[1]]->getAbundance(i);
+ if (pDataArray->groupings[(it->second)[0]][(it->second)[1]]->getAbundance(i) != 0.0) { numNotZero++; }
+ }
+
+ }
+
+ //mean abundance
+ float Aij = (totalAbund / (float) pDataArray->groupings[j].size());
+ terms.push_back(Aij);
+
+ //percentage of sites represented
+ Bij.push_back(numNotZero / (float) pDataArray->groupings[j].size());
+
+ AijDenominator += Aij;
+ }
+
+ float maxIndVal = 0.0;
+ for (int j = 0; j < terms.size(); j++) {
+ float thisAij = (terms[j] / AijDenominator); //relative abundance
+ float thisValue = thisAij * Bij[j] * 100.0;
+
+ //save largest
+ if (thisValue > maxIndVal) { maxIndVal = thisValue; }
+ }
+
+ randomIndicatorValues.push_back(maxIndVal);
+ }
+
+ ///////////////////////////////////////////////////////////////////////
+
+ for (int j = 0; j < pDataArray->indicatorValues.size(); j++) {
+ if (randomIndicatorValues[j] >= pDataArray->indicatorValues[j]) { pDataArray->pvalues[j]++; }
+ }
+ }
+
+ return 0;
+
+ }
+ catch(exception& e) {
+ pDataArray->m->errorOut(e, "IndicatorCommand", "MyIndicatorThreadFunction");
+ exit(1);
+ }
+}
+#endif
+
+
#endif
namesOfGroupCombos.push_back(groups);
}
}
-
- #if defined (__APPLE__) || (__MACH__) || (linux) || (__linux) || (__linux__) || (__unix__) || (__unix)
- if(processors == 1){
- data = driver(t, namesOfGroupCombos, 0, namesOfGroupCombos.size(), ct);
- }else{
- lines.clear();
- 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(linePair(startPos, numPairsPerProcessor));
- }
-
- data = createProcesses(t, namesOfGroupCombos, ct);
- }
- #else
- data = driver(t, namesOfGroupCombos, 0, namesOfGroupCombos.size(), ct);
- #endif
+
+ lines.clear();
+ 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(linePair(startPos, numPairsPerProcessor));
+ }
+
+ data = createProcesses(t, namesOfGroupCombos, ct);
return data;
EstOutput Parsimony::createProcesses(Tree* t, vector< vector<string> > namesOfGroupCombos, CountTable* ct) {
try {
-#if defined (__APPLE__) || (__MACH__) || (linux) || (__linux) || (__linux__) || (__unix__) || (__unix)
- int process = 1;
+ int process = 1;
vector<int> processIDS;
EstOutput results;
-
+
+#if defined (__APPLE__) || (__MACH__) || (linux) || (__linux) || (__linux__) || (__unix__) || (__unix)
+
//loop through and create all the processes you want
while (process != processors) {
int pid = fork();
in.close();
m->mothurRemove(s);
}
+#else
+ //fill in functions
+ vector<parsData*> pDataArray;
+ DWORD dwThreadIdArray[processors-1];
+ HANDLE hThreadArray[processors-1];
+ vector<CountTable*> cts;
+ vector<Tree*> trees;
+
+ //Create processor worker threads.
+ for( int i=1; i<processors; i++ ){
+ CountTable* copyCount = new CountTable();
+ copyCount->copy(ct);
+ Tree* copyTree = new Tree(copyCount);
+ copyTree->getCopy(t);
+
+ cts.push_back(copyCount);
+ trees.push_back(copyTree);
+
+ parsData* temppars = new parsData(m, lines[i].start, lines[i].num, namesOfGroupCombos, copyTree, copyCount);
+ pDataArray.push_back(temppars);
+ processIDS.push_back(i);
+
+ hThreadArray[i-1] = CreateThread(NULL, 0, MyParsimonyThreadFunction, pDataArray[i-1], 0, &dwThreadIdArray[i-1]);
+ }
+
+ results = driver(t, namesOfGroupCombos, lines[0].start, lines[0].num, ct);
+
+ //Wait until all threads have terminated.
+ WaitForMultipleObjects(processors-1, hThreadArray, TRUE, INFINITE);
+
+ //Close all thread handles and free memory allocations.
+ for(int i=0; i < pDataArray.size(); i++){
+ for (int j = 0; j < pDataArray[i]->results.size(); j++) { results.push_back(pDataArray[i]->results[j]); }
+ delete cts[i];
+ delete trees[i];
+ CloseHandle(hThreadArray[i]);
+ delete pDataArray[i];
+ }
- return results;
#endif
+ return results;
}
catch(exception& e) {
m->errorOut(e, "Parsimony", "createProcesses");
EstOutput driver(Tree*, vector< vector<string> >, int, int, CountTable*);
EstOutput createProcesses(Tree*, vector< vector<string> >, CountTable*);
};
-
/***********************************************************************/
+struct parsData {
+ int start;
+ int num;
+ MothurOut* m;
+ EstOutput results;
+ vector< vector<string> > namesOfGroupCombos;
+ Tree* t;
+ CountTable* ct;
+
+ parsData(){}
+ parsData(MothurOut* mout, int st, int en, vector< vector<string> > ngc, Tree* tree, CountTable* count) {
+ m = mout;
+ start = st;
+ num = en;
+ namesOfGroupCombos = ngc;
+ t = tree;
+ ct = count;
+ }
+};
+
+/**************************************************************************************************/
+#if defined (__APPLE__) || (__MACH__) || (linux) || (__linux) || (__linux__) || (__unix__) || (__unix)
+#else
+static DWORD WINAPI MyParsimonyThreadFunction(LPVOID lpParam){
+ parsData* pDataArray;
+ pDataArray = (parsData*)lpParam;
+ try {
+
+ pDataArray->results.resize(pDataArray->num);
+
+ Tree* copyTree = new Tree(pDataArray->ct);
+ int count = 0;
+
+ for (int h = pDataArray->start; h < (pDataArray->start+pDataArray->num); h++) {
+
+ if (pDataArray->m->control_pressed) { delete copyTree; return 0; }
+
+ int score = 0;
+
+ //groups in this combo
+ vector<string> groups = pDataArray->namesOfGroupCombos[h];
+
+ //copy users tree so that you can redo pgroups
+ copyTree->getCopy(pDataArray->t);
+
+ //create pgroups that reflect the groups the user want to use
+ for(int i=copyTree->getNumLeaves();i<copyTree->getNumNodes();i++){
+ copyTree->tree[i].pGroups = (copyTree->mergeUserGroups(i, groups));
+ }
+
+ for(int i=copyTree->getNumLeaves();i<copyTree->getNumNodes();i++){
+
+ if (pDataArray->m->control_pressed) { return 0; }
+
+ int lc = copyTree->tree[i].getLChild();
+ int rc = copyTree->tree[i].getRChild();
+
+ int iSize = copyTree->tree[i].pGroups.size();
+ int rcSize = copyTree->tree[rc].pGroups.size();
+ int lcSize = copyTree->tree[lc].pGroups.size();
+
+ //if isize are 0 then that branch is to be ignored
+ if (iSize == 0) { }
+ else if ((rcSize == 0) || (lcSize == 0)) { }
+ //if you have more groups than either of your kids then theres been a change.
+ else if(iSize > rcSize || iSize > lcSize){
+ score++;
+ }
+ }
+
+ pDataArray->results[count] = score;
+ count++;
+ }
+
+ delete copyTree;
+
+ return 0;
+
+ }
+ catch(exception& e) {
+ pDataArray->m->errorOut(e, "Parsimony", "MyParsimonyThreadFunction");
+ exit(1);
+ }
+}
+#endif
#endif
if (numSampledList.count(diversity[mGroups[j]].size()-1) == 0) { numSampledList.insert(diversity[mGroups[j]].size()-1); }
}
- #if defined (__APPLE__) || (__MACH__) || (linux) || (__linux) || (__linux__) || (__unix__) || (__unix)
- if(processors == 1){
- driver(trees[i], diversity, sumDiversity, iters, increment, randomLeaf, numSampledList, outCollect, outSum, true);
- }else{
- if (rarefy) {
- vector<int> procIters;
-
- int numItersPerProcessor = iters / processors;
-
- //divide iters between processes
- for (int h = 0; h < processors; h++) {
- if(h == processors - 1){
- numItersPerProcessor = iters - h * numItersPerProcessor;
- }
- procIters.push_back(numItersPerProcessor);
- }
-
- createProcesses(procIters, trees[i], diversity, sumDiversity, iters, increment, randomLeaf, numSampledList, outCollect, outSum);
-
- }else{ //no need to paralellize if you dont want to rarefy
- driver(trees[i], diversity, sumDiversity, iters, increment, randomLeaf, numSampledList, outCollect, outSum, true);
- }
- }
-
- #else
- driver(trees[i], diversity, sumDiversity, iters, increment, randomLeaf, numSampledList, outCollect, outSum, true);
- #endif
-
+ if (rarefy) {
+ vector<int> procIters;
+ int numItersPerProcessor = iters / processors;
+
+ //divide iters between processes
+ for (int h = 0; h < processors; h++) {
+ if(h == processors - 1){ numItersPerProcessor = iters - h * numItersPerProcessor; }
+ procIters.push_back(numItersPerProcessor);
+ }
+
+ createProcesses(procIters, trees[i], diversity, sumDiversity, iters, increment, randomLeaf, numSampledList, outCollect, outSum);
+
+ }else{ //no need to paralellize if you dont want to rarefy
+ driver(trees[i], diversity, sumDiversity, iters, increment, randomLeaf, numSampledList, outCollect, outSum, true);
+ }
+
if (rarefy) { printData(numSampledList, sumDiversity, outRare, iters); }
}
//**********************************************************************************************************************
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){
try {
- #if defined (__APPLE__) || (__MACH__) || (linux) || (__linux) || (__linux__) || (__unix__) || (__unix)
- int process = 1;
+ int process = 1;
vector<int> processIDS;
map< string, vector<float> >::iterator itSum;
-
+
+ #if defined (__APPLE__) || (__MACH__) || (linux) || (__linux) || (__linux__) || (__unix__) || (__unix)
+
//loop through and create all the processes you want
while (process != processors) {
int pid = fork();
in.close();
m->mothurRemove(inTemp);
}
+#else
+
+ //fill in functions
+ vector<phylodivData*> pDataArray;
+ DWORD dwThreadIdArray[processors-1];
+ HANDLE hThreadArray[processors-1];
+ vector<CountTable*> cts;
+ vector<Tree*> trees;
+ map<string, int> rootForGroup = getRootForGroups(t);
+
+ //Create processor worker threads.
+ for( int i=1; i<processors; i++ ){
+ CountTable* copyCount = new CountTable();
+ copyCount->copy(ct);
+ Tree* copyTree = new Tree(copyCount);
+ copyTree->getCopy(t);
+
+ cts.push_back(copyCount);
+ trees.push_back(copyTree);
+
+ map<string, vector<float> > copydiv = div;
+ map<string, vector<float> > copysumDiv = sumDiv;
+ vector<int> copyrandomLeaf = randomLeaf;
+ set<int> copynumSampledList = numSampledList;
+ map<string, int> copyRootForGrouping = rootForGroup;
+
+ phylodivData* temp = new phylodivData(m, procIters[i], copydiv, copysumDiv, copyTree, copyCount, increment, copyrandomLeaf, copynumSampledList, copyRootForGrouping);
+ pDataArray.push_back(temp);
+ processIDS.push_back(i);
+
+ hThreadArray[i-1] = CreateThread(NULL, 0, MyPhyloDivThreadFunction, pDataArray[i-1], 0, &dwThreadIdArray[i-1]);
+ }
+
+ driver(t, div, sumDiv, procIters[0], increment, randomLeaf, numSampledList, outCollect, outSum, true);
+
+ //Wait until all threads have terminated.
+ WaitForMultipleObjects(processors-1, hThreadArray, TRUE, INFINITE);
+
+ //Close all thread handles and free memory allocations.
+ for(int i=0; i < pDataArray.size(); i++){
+ for (itSum = pDataArray[i]->sumDiv.begin(); itSum != pDataArray[i]->sumDiv.end(); itSum++) {
+ for (int k = 0; k < (itSum->second).size(); k++) {
+ sumDiv[itSum->first][k] += (itSum->second)[k];
+ }
+ }
+ delete cts[i];
+ delete trees[i];
+ CloseHandle(hThreadArray[i]);
+ delete pDataArray[i];
+ }
#endif
else { score = div[mGroups[j]][numSampled] / (float)numIters; }
out << setprecision(4) << score << endl;
- cout << mGroups[j] << '\t' << numSampled << '\t'<< setprecision(4) << score << endl;
+ //cout << mGroups[j] << '\t' << numSampled << '\t'<< setprecision(4) << score << endl;
}
out.close();
};
+/***********************************************************************/
+struct phylodivData {
+ int numIters;
+ MothurOut* m;
+ map< string, vector<float> > div;
+ map<string, vector<float> > sumDiv;
+ map<string, int> rootForGroup;
+ vector<int> randomLeaf;
+ set<int> numSampledList;
+ int increment;
+ Tree* t;
+ CountTable* ct;
+ bool includeRoot;
+
+
+ phylodivData(){}
+ phylodivData(MothurOut* mout, int ni, map< string, vector<float> > cd, map< string, vector<float> > csd, Tree* tree, CountTable* count, int incre, vector<int> crl, set<int> nsl, map<string, int> rfg) {
+ m = mout;
+ t = tree;
+ ct = count;
+ div = cd;
+ numIters = ni;
+ sumDiv = csd;
+ increment = incre;
+ randomLeaf = crl;
+ numSampledList = nsl;
+ rootForGroup = rfg;
+ }
+};
+
+/**************************************************************************************************/
+#if defined (__APPLE__) || (__MACH__) || (linux) || (__linux) || (__linux__) || (__unix__) || (__unix)
+#else
+static DWORD WINAPI MyPhyloDivThreadFunction(LPVOID lpParam){
+ phylodivData* pDataArray;
+ pDataArray = (phylodivData*)lpParam;
+ try {
+ int numLeafNodes = pDataArray->randomLeaf.size();
+ vector<string> mGroups = pDataArray->m->getGroups();
+
+ for (int l = 0; l < pDataArray->numIters; l++) {
+ random_shuffle(pDataArray->randomLeaf.begin(), pDataArray->randomLeaf.end());
+
+ //initialize counts
+ map<string, int> counts;
+ vector< map<string, bool> > countedBranch;
+ for (int i = 0; i < pDataArray->t->getNumNodes(); i++) {
+ map<string, bool> temp;
+ for (int j = 0; j < mGroups.size(); j++) { temp[mGroups[j]] = false; }
+ countedBranch.push_back(temp);
+ }
+
+ for (int j = 0; j < mGroups.size(); j++) { counts[mGroups[j]] = 0; }
+
+ for(int k = 0; k < numLeafNodes; k++){
+
+ if (pDataArray->m->control_pressed) { return 0; }
+
+ //calc branch length of randomLeaf k
+ //vector<float> br = calcBranchLength(t, randomLeaf[k], countedBranch, rootForGroup);
+ //(Tree* t, int leaf, vector< map<string, bool> >& counted, map<string, int> roots
+ /////////////////////////////////////////////////////////////////////////////////////
+ vector<float> br;
+ int index = pDataArray->randomLeaf[k];
+
+ vector<string> groups = pDataArray->t->tree[pDataArray->randomLeaf[k]].getGroup();
+ br.resize(groups.size(), 0.0);
+
+ //you are a leaf
+ if(pDataArray->t->tree[index].getBranchLength() != -1){
+ for (int k = 0; k < groups.size(); k++) {
+ br[k] += abs(pDataArray->t->tree[index].getBranchLength());
+ }
+ }
+
+ index = pDataArray->t->tree[index].getParent();
+
+ //while you aren't at root
+ while(pDataArray->t->tree[index].getParent() != -1){
+
+ if (pDataArray->m->control_pressed) { return 0; }
+
+ for (int k = 0; k < groups.size(); k++) {
+
+ if (index >= pDataArray->rootForGroup[groups[k]]) { countedBranch[index][groups[k]] = true; } //if you are at this groups "root", then say we are done
+
+ if (!countedBranch[index][groups[k]]){ //if counted[index][groups[k] is true this groups has already added all br from here to root, so quit early
+ if (pDataArray->t->tree[index].getBranchLength() != -1) {
+ br[k] += abs(pDataArray->t->tree[index].getBranchLength());
+ }
+ countedBranch[index][groups[k]] = true;
+ }
+ }
+ index = pDataArray->t->tree[index].getParent();
+ }
+ /////////////////////////////////////////////////////////////////////////////////////
+
+ //for each group in the groups update the total branch length accounting for the names file
+ groups = pDataArray->t->tree[pDataArray->randomLeaf[k]].getGroup();
+
+ for (int j = 0; j < groups.size(); j++) {
+
+ if (pDataArray->m->inUsersGroups(groups[j], mGroups)) {
+ int numSeqsInGroupJ = 0;
+ map<string, int>::iterator it;
+ it = pDataArray->t->tree[pDataArray->randomLeaf[k]].pcount.find(groups[j]);
+ if (it != pDataArray->t->tree[pDataArray->randomLeaf[k]].pcount.end()) { //this leaf node contains seqs from group j
+ numSeqsInGroupJ = it->second;
+ }
+
+ if (numSeqsInGroupJ != 0) { pDataArray->div[groups[j]][(counts[groups[j]]+1)] = pDataArray->div[groups[j]][counts[groups[j]]] + br[j]; }
+
+ for (int s = (counts[groups[j]]+2); s <= (counts[groups[j]]+numSeqsInGroupJ); s++) {
+ pDataArray->div[groups[j]][s] = pDataArray->div[groups[j]][s-1]; //update counts, but don't add in redundant branch lengths
+ }
+ counts[groups[j]] += numSeqsInGroupJ;
+ }
+ }
+ }
+
+
+ //add this diversity to the sum
+ for (int j = 0; j < mGroups.size(); j++) {
+ for (int g = 0; g < pDataArray->div[mGroups[j]].size(); g++) {
+ pDataArray->sumDiv[mGroups[j]][g] += pDataArray->div[mGroups[j]][g];
+ }
+ }
+
+ }
+
+ return 0;
+ }
+ catch(exception& e) {
+ pDataArray->m->errorOut(e, "PhyloDiversityCommand", "MyPhyloDivThreadFunction");
+ exit(1);
+ }
+}
+#endif
+
#endif
while(fileHandle && lt == 1){ //let's assume it's a triangular matrix...
- fileHandle >> firstName >> secondName >> distance; // get the row and column names and distance
+ fileHandle >> firstName; m->gobble(fileHandle);
+ fileHandle >> secondName; m->gobble(fileHandle);
+ fileHandle >> distance; // get the row and column names and distance
+
+ if (m->debug) { cout << firstName << '\t' << secondName << '\t' << distance << endl; }
if (m->control_pressed) { fileHandle.close(); delete reading; return 0; }
m->openInputFile(distFile, fileHandle); //let's start over
while(fileHandle){
- fileHandle >> firstName >> secondName >> distance;
+ fileHandle >> firstName; m->gobble(fileHandle);
+ fileHandle >> secondName; m->gobble(fileHandle);
+ fileHandle >> distance; // get the row and column names and distance
if (m->control_pressed) { fileHandle.close(); delete reading; return 0; }
while(fileHandle && lt == 1){ //let's assume it's a triangular matrix...
- fileHandle >> firstName >> secondName >> distance; // get the row and column names and distance
-
+ fileHandle >> firstName; m->gobble(fileHandle);
+ fileHandle >> secondName; m->gobble(fileHandle);
+ fileHandle >> distance; // get the row and column names and distance
+
if (m->control_pressed) { fileHandle.close(); delete reading; return 0; }
int itA = countTable->get(firstName);
m->openInputFile(distFile, fileHandle); //let's start over
while(fileHandle){
- fileHandle >> firstName >> secondName >> distance;
+ fileHandle >> firstName; m->gobble(fileHandle);
+ fileHandle >> secondName; m->gobble(fileHandle);
+ fileHandle >> distance; // get the row and column names and distance
if (m->control_pressed) { fileHandle.close(); delete reading; return 0; }
int numGroups = 0;
- ofstream outFile;
+ //ofstream outFile;
ifstream dFile;
m->openInputFile(distFile, dFile);
//have we reached the max buffer size
if (numOutputs[groupID] > 60) { //write out sequence
+ ofstream outFile;
outFile.open(fileName.c_str(), ios::app);
outFile << outputs[groupID] << seqA << '\t' << seqB << '\t' << dist << endl;
outFile.close();
//if groupB is written to file it is above buffer size so read and write to new merged file
if (wroteOutPut[groupIDB]) {
string fileName2 = distFile + "." + toString(groupIDB) + ".temp";
- ifstream fileB(fileName2.c_str(), ios::ate);
+ /*ifstream fileB(fileName2.c_str(), ios::ate);
outFile.open(fileName.c_str(), ios::app);
outFile << temp.substr(0, lastRead);
delete memblock;
- fileB.close();
+ fileB.close();*/
+ m->appendFiles(fileName2, fileName);
m->mothurRemove(fileName2);
+
//write out the merged memory
if (numOutputs[groupID] > 60) {
- outFile << outputs[groupID];
+ ofstream tempOut;
+ m->openOutputFile(fileName, tempOut);
+ tempOut << outputs[groupID];
outputs[groupID] = "";
numOutputs[groupID] = 0;
+ tempOut.close();
}
- outFile.close();
+ //outFile.close();
wroteOutPut[groupID] = true;
wroteOutPut[groupIDB] = false;
if (wroteOutPut[groupIDA]) {
string fileName2 = distFile + "." + toString(groupIDA) + ".temp";
- ifstream fileB(fileName2.c_str(), ios::ate);
+ /*ifstream fileB(fileName2.c_str(), ios::ate);
outFile.open(fileName.c_str(), ios::app);
delete memblock;
- fileB.close();
+ fileB.close();*/
+ m->appendFiles(fileName2, fileName);
m->mothurRemove(fileName2);
//write out the merged memory
if (numOutputs[groupID] > 60) {
- outFile << outputs[groupID];
+ ofstream tempOut;
+ m->openOutputFile(fileName, tempOut);
+ tempOut << outputs[groupID];
outputs[groupID] = "";
numOutputs[groupID] = 0;
+ tempOut.close();
}
- outFile.close();
+ //outFile.close();
wroteOutPut[groupID] = true;
wroteOutPut[groupIDA] = false;
//remove old names files just in case
if (numOutputs[i] > 0) {
+ ofstream outFile;
outFile.open(fileName.c_str(), ios::app);
outFile << outputs[i];
outFile.close();
vector<summarySharedData*> pDataArray;
DWORD dwThreadIdArray[processors-1];
- HANDLE hThreadArray[processors-1];
+ HANDLE hThreadArray[processors-1];
//Create processor worker threads.
for( int i=1; i<processors; i++ ){
temp->setGroup(thisLookup[k]->getGroup());
newLookup.push_back(temp);
}
+
//for each bin
for (int k = 0; k < thisLookup[0]->getNumBins(); k++) {
namesOfGroupCombos.push_back(groups);
}
}
+
+ lines.clear();
+ 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(linePair(startPos, numPairsPerProcessor));
+ }
+
+ data = createProcesses(t, namesOfGroupCombos, ct);
- #if defined (__APPLE__) || (__MACH__) || (linux) || (__linux) || (__linux__) || (__unix__) || (__unix)
- if(processors == 1){
- data = driver(t, namesOfGroupCombos, 0, namesOfGroupCombos.size(), ct);
- }else{
- 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(linePair(startPos, numPairsPerProcessor));
- }
- data = createProcesses(t, namesOfGroupCombos, ct);
- lines.clear();
- }
- #else
- data = driver(t, namesOfGroupCombos, 0, namesOfGroupCombos.size(), ct);
- #endif
-
+ lines.clear();
+
return data;
}
catch(exception& e) {
EstOutput Unweighted::createProcesses(Tree* t, vector< vector<string> > namesOfGroupCombos, CountTable* ct) {
try {
-#if defined (__APPLE__) || (__MACH__) || (linux) || (__linux) || (__linux__) || (__unix__) || (__unix)
- int process = 1;
+ int process = 1;
vector<int> processIDS;
EstOutput results;
+#if defined (__APPLE__) || (__MACH__) || (linux) || (__linux) || (__linux__) || (__unix__) || (__unix)
+
//loop through and create all the processes you want
while (process != processors) {
in.close();
m->mothurRemove(s);
}
+#else
+ //fill in functions
+ vector<unweightedData*> pDataArray;
+ DWORD dwThreadIdArray[processors-1];
+ HANDLE hThreadArray[processors-1];
+ vector<CountTable*> cts;
+ vector<Tree*> trees;
+
+ //Create processor worker threads.
+ for( int i=1; i<processors; i++ ){
+ CountTable* copyCount = new CountTable();
+ copyCount->copy(ct);
+ Tree* copyTree = new Tree(copyCount);
+ copyTree->getCopy(t);
+
+ cts.push_back(copyCount);
+ trees.push_back(copyTree);
+
+ unweightedData* tempweighted = new unweightedData(m, lines[i].start, lines[i].num, namesOfGroupCombos, copyTree, copyCount, includeRoot);
+ pDataArray.push_back(tempweighted);
+ processIDS.push_back(i);
+
+ hThreadArray[i-1] = CreateThread(NULL, 0, MyUnWeightedThreadFunction, pDataArray[i-1], 0, &dwThreadIdArray[i-1]);
+ }
+
+ results = driver(t, namesOfGroupCombos, lines[0].start, lines[0].num, ct);
- //m->mothurOut("DONE."); m->mothurOutEndLine(); m->mothurOutEndLine();
+ //Wait until all threads have terminated.
+ WaitForMultipleObjects(processors-1, hThreadArray, TRUE, INFINITE);
- return results;
-#endif
+ //Close all thread handles and free memory allocations.
+ for(int i=0; i < pDataArray.size(); i++){
+ for (int j = 0; j < pDataArray[i]->results.size(); j++) { results.push_back(pDataArray[i]->results[j]); }
+ delete cts[i];
+ delete trees[i];
+ CloseHandle(hThreadArray[i]);
+ delete pDataArray[i];
+ }
+
+#endif
+ return results;
}
catch(exception& e) {
m->errorOut(e, "Unweighted", "createProcesses");
int count = 0;
int total = num;
- int twentyPercent = (total * 0.20);
- if (twentyPercent == 0) { twentyPercent = 1; }
-
-
+
for (int h = start; h < (start+num); h++) {
if (m->control_pressed) { return results; }
}
count++;
- //report progress
- //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(); }
- }
-
- //report progress
- //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(); }
+ }
return results;
}
namesOfGroupCombos.push_back(groups);
}
}
+
+ lines.clear();
+ 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(linePair(startPos, numPairsPerProcessor));
+ }
+
+ data = createProcesses(t, namesOfGroupCombos, true, ct);
+ lines.clear();
- #if defined (__APPLE__) || (__MACH__) || (linux) || (__linux) || (__linux__) || (__unix__) || (__unix)
- if(processors == 1){
- data = driver(t, namesOfGroupCombos, 0, namesOfGroupCombos.size(), true, ct);
- }else{
- 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(linePair(startPos, numPairsPerProcessor));
- }
-
- data = createProcesses(t, namesOfGroupCombos, true, ct);
-
- lines.clear();
- }
- #else
- data = driver(t, namesOfGroupCombos, 0, namesOfGroupCombos.size(), true, ct);
- #endif
-
return data;
}
catch(exception& e) {
EstOutput Unweighted::createProcesses(Tree* t, vector< vector<string> > namesOfGroupCombos, bool usingGroups, CountTable* ct) {
try {
-#if defined (__APPLE__) || (__MACH__) || (linux) || (__linux) || (__linux__) || (__unix__) || (__unix)
- int process = 1;
+ int process = 1;
vector<int> processIDS;
EstOutput results;
-
+#if defined (__APPLE__) || (__MACH__) || (linux) || (__linux) || (__linux__) || (__unix__) || (__unix)
+
//loop through and create all the processes you want
while (process != processors) {
int pid = fork();
in.close();
m->mothurRemove(s);
}
+#else
+ //for some reason it doesn't seem to be calculating hte random trees scores. all scores are the same even though copytree appears to be randomized.
+
+ /*
+ //fill in functions
+ vector<unweightedData*> pDataArray;
+ DWORD dwThreadIdArray[processors-1];
+ HANDLE hThreadArray[processors-1];
+ vector<CountTable*> cts;
+ vector<Tree*> trees;
+
+ //Create processor worker threads.
+ for( int i=1; i<processors; i++ ){
+ CountTable* copyCount = new CountTable();
+ copyCount->copy(ct);
+ Tree* copyTree = new Tree(copyCount);
+ copyTree->getCopy(t);
+
+ cts.push_back(copyCount);
+ trees.push_back(copyTree);
+
+ unweightedData* tempweighted = new unweightedData(m, lines[i].start, lines[i].num, namesOfGroupCombos, copyTree, copyCount, includeRoot);
+ pDataArray.push_back(tempweighted);
+ processIDS.push_back(i);
+
+ hThreadArray[i-1] = CreateThread(NULL, 0, MyUnWeightedRandomThreadFunction, pDataArray[i-1], 0, &dwThreadIdArray[i-1]);
+ }
+
+ results = driver(t, namesOfGroupCombos, lines[0].start, lines[0].num, usingGroups, ct);
- return results;
-#endif
+ //Wait until all threads have terminated.
+ WaitForMultipleObjects(processors-1, hThreadArray, TRUE, INFINITE);
+
+ //Close all thread handles and free memory allocations.
+ for(int i=0; i < pDataArray.size(); i++){
+ for (int j = 0; j < pDataArray[i]->results.size(); j++) { results.push_back(pDataArray[i]->results[j]); }
+ delete cts[i];
+ delete trees[i];
+ CloseHandle(hThreadArray[i]);
+ delete pDataArray[i];
+ } */
+
+ results = driver(t, namesOfGroupCombos, 0, namesOfGroupCombos.size(), usingGroups, ct);
+#endif
+ return results;
}
catch(exception& e) {
m->errorOut(e, "Unweighted", "createProcesses");
if (isnan(UW) || isinf(UW)) { UW = 0; }
results[count] = UW;
- }
+ }
count++;
}
};
/***********************************************************************/
+struct unweightedData {
+ int start;
+ int num;
+ MothurOut* m;
+ EstOutput results;
+ vector< vector<string> > namesOfGroupCombos;
+ Tree* t;
+ CountTable* ct;
+ bool includeRoot;
+
+ unweightedData(){}
+ unweightedData(MothurOut* mout, int st, int en, vector< vector<string> > ngc, Tree* tree, CountTable* count, bool ir) {
+ m = mout;
+ start = st;
+ num = en;
+ namesOfGroupCombos = ngc;
+ t = tree;
+ ct = count;
+ includeRoot = ir;
+ }
+};
+
+/**************************************************************************************************/
+#if defined (__APPLE__) || (__MACH__) || (linux) || (__linux) || (__linux__) || (__unix__) || (__unix)
+#else
+static DWORD WINAPI MyUnWeightedThreadFunction(LPVOID lpParam){
+ unweightedData* pDataArray;
+ pDataArray = (unweightedData*)lpParam;
+ try {
+ pDataArray->results.resize(pDataArray->num);
+ map< vector<string>, set<int> > rootForGrouping;
+
+ int count = 0;
+
+ for (int h = pDataArray->start; h < (pDataArray->start+pDataArray->num); h++) {
+ if (pDataArray->m->control_pressed) { return 0; }
+
+ double UniqueBL=0.0000; //a branch length is unique if it's chidren are from the same group
+ double totalBL = 0.00; //all branch lengths
+ double UW = 0.00; //Unweighted Value = UniqueBL / totalBL;
+
+ //find a node that belongs to one of the groups in this combo
+ int nodeBelonging = -1;
+ for (int g = 0; g < pDataArray->namesOfGroupCombos[h].size(); g++) {
+ if (pDataArray->t->groupNodeInfo[pDataArray->namesOfGroupCombos[h][g]].size() != 0) { nodeBelonging = pDataArray->t->groupNodeInfo[pDataArray->namesOfGroupCombos[h][g]][0]; break; }
+ }
+
+ //sanity check
+ if (nodeBelonging == -1) {
+ pDataArray->m->mothurOut("[WARNING]: cannot find a nodes in the tree from grouping ");
+ for (int g = 0; g < pDataArray->namesOfGroupCombos[h].size()-1; g++) { pDataArray->m->mothurOut(pDataArray->namesOfGroupCombos[h][g] + "-"); }
+ pDataArray->m->mothurOut(pDataArray->namesOfGroupCombos[h][pDataArray->namesOfGroupCombos[h].size()-1]);
+ pDataArray->m->mothurOut(", skipping."); pDataArray->m->mothurOutEndLine(); pDataArray->results[count] = UW;
+ }else{
+
+ //if including the root this clears rootForGrouping[namesOfGroupCombos[h]]
+ //getRoot(t, nodeBelonging, namesOfGroupCombos[h]);
+ /////////////////////////////////////////////////////////////////////////////
+ //you are a leaf so get your parent
+ vector<string> grouping = pDataArray->namesOfGroupCombos[h];
+ int index = pDataArray->t->tree[nodeBelonging].getParent();
+
+ if (pDataArray->includeRoot) {
+ rootForGrouping[grouping].clear();
+ }else {
+
+ //my parent is a potential root
+ rootForGrouping[grouping].insert(index);
+
+ //while you aren't at root
+ while(pDataArray->t->tree[index].getParent() != -1){
+ //cout << index << endl;
+ if (pDataArray->m->control_pressed) { return 0; }
+
+ //am I the root for this grouping? if so I want to stop "early"
+ //does my sibling have descendants from the users groups?
+ //if so I am not the root
+ int parent = pDataArray->t->tree[index].getParent();
+ int lc = pDataArray->t->tree[parent].getLChild();
+ int rc = pDataArray->t->tree[parent].getRChild();
+
+ int sib = lc;
+ if (lc == index) { sib = rc; }
+
+ map<string, int>::iterator itGroup;
+ int pcountSize = 0;
+ for (int j = 0; j < grouping.size(); j++) {
+ map<string, int>::iterator itGroup = pDataArray->t->tree[sib].pcount.find(grouping[j]);
+ if (itGroup != pDataArray->t->tree[sib].pcount.end()) { pcountSize++; if (pcountSize > 1) { break; } }
+ }
+
+ //if yes, I am not the root
+ if (pcountSize != 0) {
+ rootForGrouping[grouping].clear();
+ rootForGrouping[grouping].insert(parent);
+ }
+
+ index = parent;
+ }
+
+ //get all nodes above the root to add so we don't add their u values above
+ index = *(rootForGrouping[grouping].begin());
+ while(pDataArray->t->tree[index].getParent() != -1){
+ int parent = pDataArray->t->tree[index].getParent();
+ rootForGrouping[grouping].insert(parent);
+ //cout << parent << " in root" << endl;
+ index = parent;
+ }
+ }
+ /////////////////////////////////////////////////////////////////////////////
+
+ for(int i=0;i<pDataArray->t->getNumNodes();i++){
+
+ if (pDataArray->m->control_pressed) { return 0; }
+ //cout << i << endl;
+ //pcountSize = 0, they are from a branch that is entirely from a group the user doesn't want
+ //pcountSize = 2, not unique to one group
+ //pcountSize = 1, unique to one group
+
+ int pcountSize = 0;
+ for (int j = 0; j < pDataArray->namesOfGroupCombos[h].size(); j++) {
+ map<string, int>::iterator itGroup = pDataArray->t->tree[i].pcount.find(pDataArray->namesOfGroupCombos[h][j]);
+ if (itGroup != pDataArray->t->tree[i].pcount.end()) { pcountSize++; if (pcountSize > 1) { break; } }
+ }
+
+
+ //unique calc
+ if (pcountSize == 0) { }
+ else if ((pDataArray->t->tree[i].getBranchLength() != -1) && (pcountSize == 1) && (rootForGrouping[pDataArray->namesOfGroupCombos[h]].count(i) == 0)) { //you have a unique branch length and you are not the root
+ UniqueBL += abs(pDataArray->t->tree[i].getBranchLength());
+ }
+
+ //total calc
+ if (pcountSize == 0) { }
+ else if ((pDataArray->t->tree[i].getBranchLength() != -1) && (pcountSize != 0) && (rootForGrouping[pDataArray->namesOfGroupCombos[h]].count(i) == 0)) { //you have a branch length and you are not the root
+ totalBL += abs(pDataArray->t->tree[i].getBranchLength());
+ }
+ }
+ //cout << UniqueBL << '\t' << totalBL << endl;
+ UW = (UniqueBL / totalBL);
+
+ if (isnan(UW) || isinf(UW)) { UW = 0; }
+
+ pDataArray->results[count] = UW;
+ }
+ count++;
+ }
+
+ return 0;
+
+ }
+ catch(exception& e) {
+ pDataArray->m->errorOut(e, "UnWeighted", "MyUnWeightedThreadFunction");
+ exit(1);
+ }
+}
+/**************************************************************************************************/
+
+static DWORD WINAPI MyUnWeightedRandomThreadFunction(LPVOID lpParam){
+ unweightedData* pDataArray;
+ pDataArray = (unweightedData*)lpParam;
+ try {
+ pDataArray->results.resize(pDataArray->num);
+
+ int count = 0;
+
+ Tree* copyTree = new Tree(pDataArray->ct);
+
+ for (int h = pDataArray->start; h < (pDataArray->start+pDataArray->num); h++) {
+
+ if (pDataArray->m->control_pressed) { return 0; }
+
+ map< vector<string>, set<int> > rootForGrouping;
+
+ //copy random tree passed in
+ copyTree->getCopy(pDataArray->t);
+
+ //swap labels in the groups you want to compare
+ copyTree->assembleRandomUnifracTree(pDataArray->namesOfGroupCombos[h]);
+
+ double UniqueBL=0.0000; //a branch length is unique if it's chidren are from the same group
+ double totalBL = 0.00; //all branch lengths
+ double UW = 0.00; //Unweighted Value = UniqueBL / totalBL;
+ //find a node that belongs to one of the groups in this combo
+ int nodeBelonging = -1;
+ for (int g = 0; g < pDataArray->namesOfGroupCombos[h].size(); g++) {
+ if (copyTree->groupNodeInfo[pDataArray->namesOfGroupCombos[h][g]].size() != 0) { nodeBelonging = copyTree->groupNodeInfo[pDataArray->namesOfGroupCombos[h][g]][0]; break; }
+ }
+
+ //sanity check
+ if (nodeBelonging == -1) {
+ pDataArray->m->mothurOut("[WARNING]: cannot find a nodes in the tree from grouping ");
+ for (int g = 0; g < pDataArray->namesOfGroupCombos[h].size()-1; g++) { pDataArray->m->mothurOut(pDataArray->namesOfGroupCombos[h][g] + "-"); }
+ pDataArray->m->mothurOut(pDataArray->namesOfGroupCombos[h][pDataArray->namesOfGroupCombos[h].size()-1]);
+ pDataArray->m->mothurOut(", skipping."); pDataArray->m->mothurOutEndLine(); pDataArray->results[count] = UW;
+ }else{
+
+ //if including the root this clears rootForGrouping[namesOfGroupCombos[h]]
+ //getRoot(copyTree, nodeBelonging, namesOfGroupCombos[h]);
+ /////////////////////////////////////////////////////////////////////////////
+ //you are a leaf so get your parent
+ vector<string> grouping = pDataArray->namesOfGroupCombos[h];
+ int index = copyTree->tree[nodeBelonging].getParent();
+
+ if (pDataArray->includeRoot) {
+ rootForGrouping[grouping].clear();
+ }else {
+
+ //my parent is a potential root
+ rootForGrouping[grouping].insert(index);
+
+ //while you aren't at root
+ while(copyTree->tree[index].getParent() != -1){
+ //cout << index << endl;
+ if (pDataArray->m->control_pressed) { return 0; }
+
+ //am I the root for this grouping? if so I want to stop "early"
+ //does my sibling have descendants from the users groups?
+ //if so I am not the root
+ int parent = copyTree->tree[index].getParent();
+ int lc = copyTree->tree[parent].getLChild();
+ int rc = copyTree->tree[parent].getRChild();
+
+ int sib = lc;
+ if (lc == index) { sib = rc; }
+
+ map<string, int>::iterator itGroup;
+ int pcountSize = 0;
+ for (int j = 0; j < grouping.size(); j++) {
+ map<string, int>::iterator itGroup = copyTree->tree[sib].pcount.find(grouping[j]);
+ if (itGroup != copyTree->tree[sib].pcount.end()) { pcountSize++; if (pcountSize > 1) { break; } }
+ }
+
+ //if yes, I am not the root
+ if (pcountSize != 0) {
+ rootForGrouping[grouping].clear();
+ rootForGrouping[grouping].insert(parent);
+ }
+
+ index = parent;
+ }
+
+ //get all nodes above the root to add so we don't add their u values above
+ index = *(rootForGrouping[grouping].begin());
+ while(copyTree->tree[index].getParent() != -1){
+ int parent = copyTree->tree[index].getParent();
+ rootForGrouping[grouping].insert(parent);
+ //cout << parent << " in root" << endl;
+ index = parent;
+ }
+ }
+ /////////////////////////////////////////////////////////////////////////////
+ for(int i=0;i<copyTree->getNumNodes();i++){
+
+ if (pDataArray->m->control_pressed) { return 0; }
+
+ //pcountSize = 0, they are from a branch that is entirely from a group the user doesn't want
+ //pcountSize = 2, not unique to one group
+ //pcountSize = 1, unique to one group
+ int pcountSize = 0;
+ for (int j = 0; j < pDataArray->namesOfGroupCombos[h].size(); j++) {
+ map<string, int>::iterator itGroup = copyTree->tree[i].pcount.find(pDataArray->namesOfGroupCombos[h][j]);
+ if (itGroup != copyTree->tree[i].pcount.end()) { pcountSize++; if (pcountSize > 1) { break; } }
+ }
+
+ //unique calc
+ if (pcountSize == 0) { }
+ else if ((copyTree->tree[i].getBranchLength() != -1) && (pcountSize == 1) && (rootForGrouping[pDataArray->namesOfGroupCombos[h]].count(i) == 0)) { //you have a unique branch length and you are not the root
+ UniqueBL += abs(copyTree->tree[i].getBranchLength());
+ }
+
+ //total calc
+ if (pcountSize == 0) { }
+ else if ((copyTree->tree[i].getBranchLength() != -1) && (pcountSize != 0) && (rootForGrouping[pDataArray->namesOfGroupCombos[h]].count(i) == 0)) { //you have a branch length and you are not the root
+ totalBL += abs(copyTree->tree[i].getBranchLength());
+ }
+
+ }
+ cout << h << '\t' << UniqueBL << '\t' << totalBL << endl;
+ UW = (UniqueBL / totalBL);
+
+ if (isnan(UW) || isinf(UW)) { UW = 0; }
+
+ pDataArray->results[count] = UW;
+ cout << h << '\t' << UW << endl;
+ }
+ count++;
+
+ }
+
+ delete copyTree;
+
+ return 0;
+ }
+ catch(exception& e) {
+ pDataArray->m->errorOut(e, "UnWeighted", "MyUnWeightedRandomThreadFunction");
+ exit(1);
+ }
+}
+
+#endif
+
#endif
projects / mothur.git / commitdiff