X-Git-Url: https://git.donarmstrong.com/?p=mothur.git;a=blobdiff_plain;f=unifracweightedcommand.cpp;h=3b0c53b7437111cfe5045663a8c78370274b939e;hp=cbec7490bfe028ea02de0b4a1e024b1ce2d89e1f;hb=050a3ff02473a3d4c0980964e1a9ebe52e55d6b8;hpb=6c2b1e530a5c0bb87040e58a3e410097acdfcc3d diff --git a/unifracweightedcommand.cpp b/unifracweightedcommand.cpp index cbec749..3b0c53b 100644 --- a/unifracweightedcommand.cpp +++ b/unifracweightedcommand.cpp @@ -15,20 +15,20 @@ //********************************************************************************************************************** vector UnifracWeightedCommand::setParameters(){ try { - CommandParameter ptree("tree", "InputTypes", "", "", "none", "none", "none",false,true); parameters.push_back(ptree); - CommandParameter pname("name", "InputTypes", "", "", "NameCount", "none", "none",false,false); parameters.push_back(pname); - CommandParameter pcount("count", "InputTypes", "", "", "NameCount-CountGroup", "none", "none",false,false); parameters.push_back(pcount); - CommandParameter pgroup("group", "InputTypes", "", "", "CountGroup", "none", "none",false,false); parameters.push_back(pgroup); - CommandParameter pgroups("groups", "String", "", "", "", "", "",false,false); parameters.push_back(pgroups); - CommandParameter piters("iters", "Number", "", "1000", "", "", "",false,false); parameters.push_back(piters); - CommandParameter pprocessors("processors", "Number", "", "1", "", "", "",false,false); parameters.push_back(pprocessors); - CommandParameter psubsample("subsample", "String", "", "", "", "", "",false,false); parameters.push_back(psubsample); - CommandParameter pconsensus("consensus", "Boolean", "", "F", "", "", "",false,false); parameters.push_back(pconsensus); - CommandParameter prandom("random", "Boolean", "", "F", "", "", "",false,false); parameters.push_back(prandom); - CommandParameter pdistance("distance", "Multiple", "column-lt-square-phylip", "column", "", "", "",false,false); parameters.push_back(pdistance); - CommandParameter proot("root", "Boolean", "F", "", "", "", "",false,false); parameters.push_back(proot); - CommandParameter pinputdir("inputdir", "String", "", "", "", "", "",false,false); parameters.push_back(pinputdir); - CommandParameter poutputdir("outputdir", "String", "", "", "", "", "",false,false); parameters.push_back(poutputdir); + CommandParameter ptree("tree", "InputTypes", "", "", "none", "none", "none","weighted-wsummary",false,true,true); parameters.push_back(ptree); + CommandParameter pname("name", "InputTypes", "", "", "NameCount", "none", "none","",false,false,true); parameters.push_back(pname); + CommandParameter pcount("count", "InputTypes", "", "", "NameCount-CountGroup", "none", "none","",false,false,true); parameters.push_back(pcount); + CommandParameter pgroup("group", "InputTypes", "", "", "CountGroup", "none", "none","",false,false,true); parameters.push_back(pgroup); + CommandParameter pgroups("groups", "String", "", "", "", "", "","",false,false); parameters.push_back(pgroups); + CommandParameter piters("iters", "Number", "", "1000", "", "", "","",false,false); parameters.push_back(piters); + CommandParameter pprocessors("processors", "Number", "", "1", "", "", "","",false,false,true); parameters.push_back(pprocessors); + CommandParameter psubsample("subsample", "String", "", "", "", "", "","",false,false); parameters.push_back(psubsample); + CommandParameter pconsensus("consensus", "Boolean", "", "F", "", "", "","tree",false,false); parameters.push_back(pconsensus); + CommandParameter prandom("random", "Boolean", "", "F", "", "", "","",false,false); parameters.push_back(prandom); + CommandParameter pdistance("distance", "Multiple", "column-lt-square-phylip", "column", "", "", "","phylip-column",false,false); parameters.push_back(pdistance); + CommandParameter proot("root", "Boolean", "F", "", "", "", "","",false,false); parameters.push_back(proot); + CommandParameter pinputdir("inputdir", "String", "", "", "", "", "","",false,false); parameters.push_back(pinputdir); + CommandParameter poutputdir("outputdir", "String", "", "", "", "", "","",false,false); parameters.push_back(poutputdir); vector myArray; for (int i = 0; i < parameters.size(); i++) { myArray.push_back(parameters[i].name); } @@ -65,28 +65,22 @@ string UnifracWeightedCommand::getHelpString(){ } } //********************************************************************************************************************** -string UnifracWeightedCommand::getOutputFileNameTag(string type, string inputName=""){ - try { - string outputFileName = ""; - map >::iterator it; +string UnifracWeightedCommand::getOutputPattern(string type) { + try { + string pattern = ""; + if (type == "weighted") { pattern = "[filename],weighted-[filename],[tag],weighted"; } + else if (type == "wsummary") { pattern = "[filename],wsummary"; } + else if (type == "phylip") { pattern = "[filename],[tag],[tag2],dist"; } + else if (type == "column") { pattern = "[filename],[tag],[tag2],dist"; } + else if (type == "tree") { pattern = "[filename],[tag],[tag2],tre"; } + else { m->mothurOut("[ERROR]: No definition for type " + type + " output pattern.\n"); m->control_pressed = true; } - //is this a type this command creates - it = outputTypes.find(type); - if (it == outputTypes.end()) { m->mothurOut("[ERROR]: this command doesn't create a " + type + " output file.\n"); } - else { - if (type == "weighted") { outputFileName = "weighted"; } - else if (type == "wsummary") { outputFileName = "wsummary"; } - else if (type == "phylip") { outputFileName = "dist"; } - else if (type == "column") { outputFileName = "dist"; } - else if (type == "tree") { outputFileName = "tre"; } - else { m->mothurOut("[ERROR]: No definition for type " + type + " output file tag.\n"); m->control_pressed = true; } - } - return outputFileName; - } - catch(exception& e) { - m->errorOut(e, "UnifracWeightedCommand", "getOutputFileNameTag"); - exit(1); - } + return pattern; + } + catch(exception& e) { + m->errorOut(e, "UnifracWeightedCommand", "getOutputPattern"); + exit(1); + } } //********************************************************************************************************************** UnifracWeightedCommand::UnifracWeightedCommand(){ @@ -296,8 +290,10 @@ int UnifracWeightedCommand::execute() { delete reader; if (m->control_pressed) { delete ct; for (int i = 0; i < T.size(); i++) { delete T[i]; } return 0; } - - sumFile = outputDir + m->getSimpleName(treefile) + getOutputFileNameTag("wsummary"); + + map variables; + variables["[filename]"] = outputDir + m->getSimpleName(treefile); + sumFile = getOutputFileName("wsummary",variables); m->openOutputFile(sumFile, outSum); outputNames.push_back(sumFile); outputTypes["wsummary"].push_back(sumFile); @@ -305,7 +301,8 @@ int UnifracWeightedCommand::execute() { string s; //to make work with setgroups Groups = m->getGroups(); vector nameGroups = ct->getNamesOfGroups(); - util.setGroups(Groups, nameGroups, s, numGroups, "weighted"); //sets the groups the user wants to analyze + if (nameGroups.size() < 2) { m->mothurOut("[ERROR]: You cannot run unifrac.weighted with less than 2 groups, aborting.\n"); delete ct; for (int i = 0; i < T.size(); i++) { delete T[i]; } return 0; } + util.setGroups(Groups, nameGroups, s, numGroups, "weighted"); //sets the groups the user wants to analyze m->setGroups(Groups); if (m->control_pressed) { delete ct; for (int i = 0; i < T.size(); i++) { delete T[i]; } return 0; } @@ -358,9 +355,11 @@ int UnifracWeightedCommand::execute() { vector randomData; randomData.resize(numComp,0); //weighted score info for random trees. data[0] = weightedscore AB, data[1] = weightedscore AC... if (random) { - output = new ColumnFile(outputDir + m->getSimpleName(treefile) + toString(i+1) + "." + getOutputFileNameTag("weighted"), itersString); - outputNames.push_back(outputDir + m->getSimpleName(treefile) + toString(i+1) + "." + getOutputFileNameTag("weighted")); - outputTypes["weighted"].push_back(outputDir + m->getSimpleName(treefile) + toString(i+1) + "." + getOutputFileNameTag("weighted")); + variables["[filename]"] = outputDir + m->getSimpleName(treefile); + variables["[tag]"] = toString(i+1); + string wFileName = getOutputFileName("weighted", variables); + output = new ColumnFile(wFileName, itersString); + outputNames.push_back(wFileName); outputTypes["wweighted"].push_back(wFileName); } userData = weighted.getValues(T[i], processors, outputDir); //userData[0] = weightedscore @@ -384,16 +383,20 @@ int UnifracWeightedCommand::execute() { //subsample loop vector< vector > calcDistsTotals; //each iter, each groupCombos dists. this will be used to make .dist files for (int thisIter = 0; thisIter < subsampleIters; thisIter++) { //subsampleIters=0, if subsample=f. - if (m->control_pressed) { break; } //copy to preserve old one - would do this in subsample but memory cleanup becomes messy. CountTable* newCt = new CountTable(); + int sampleTime = 0; + if (m->debug) { sampleTime = time(NULL); } + //uses method of setting groups to doNotIncludeMe SubSample sample; Tree* subSampleTree = sample.getSample(T[i], ct, newCt, subsampleSize); - + + if (m->debug) { m->mothurOut("[DEBUG]: iter " + toString(thisIter) + " took " + toString(time(NULL) - sampleTime) + " seconds to sample tree.\n"); } + //call new weighted function vector iterData; iterData.resize(numComp,0); Weighted thisWeighted(includeRoot); @@ -405,7 +408,7 @@ int UnifracWeightedCommand::execute() { delete newCt; delete subSampleTree; - if((thisIter+1) % 100 == 0){ m->mothurOut(toString(thisIter+1)); m->mothurOutEndLine(); } + if((thisIter+1) % 100 == 0){ m->mothurOutJustToScreen(toString(thisIter+1)+"\n"); } } if (m->control_pressed) { delete ct; for (int i = 0; i < T.size(); i++) { delete T[i]; } if (random) { delete output; } outSum.close(); for (int i = 0; i < outputNames.size(); i++) { m->mothurRemove(outputNames[i]); } return 0; } @@ -462,40 +465,28 @@ int UnifracWeightedCommand::getAverageSTDMatrices(vector< vector >& dist //we need to find the average distance and standard deviation for each groups distance //finds sum - vector averages; averages.resize(numComp, 0); - for (int thisIter = 0; thisIter < subsampleIters; thisIter++) { - for (int i = 0; i < dists[thisIter].size(); i++) { - averages[i] += dists[thisIter][i]; - } - } - - //finds average. - for (int i = 0; i < averages.size(); i++) { averages[i] /= (float) subsampleIters; } + vector averages = m->getAverages(dists); //find standard deviation - vector stdDev; stdDev.resize(numComp, 0); - - for (int thisIter = 0; thisIter < iters; thisIter++) { //compute the difference of each dist from the mean, and square the result of each - for (int j = 0; j < dists[thisIter].size(); j++) { - stdDev[j] += ((dists[thisIter][j] - averages[j]) * (dists[thisIter][j] - averages[j])); - } - } - for (int i = 0; i < stdDev.size(); i++) { - stdDev[i] /= (float) subsampleIters; - stdDev[i] = sqrt(stdDev[i]); - } + vector stdDev = m->getStandardDeviation(dists, averages); //make matrix with scores in it - vector< vector > avedists; avedists.resize(m->getNumGroups()); + vector< vector > avedists; //avedists.resize(m->getNumGroups()); for (int i = 0; i < m->getNumGroups(); i++) { - avedists[i].resize(m->getNumGroups(), 0.0); + vector temp; + for (int j = 0; j < m->getNumGroups(); j++) { temp.push_back(0.0); } + avedists.push_back(temp); } //make matrix with scores in it - vector< vector > stddists; stddists.resize(m->getNumGroups()); + vector< vector > stddists; //stddists.resize(m->getNumGroups()); for (int i = 0; i < m->getNumGroups(); i++) { - stddists[i].resize(m->getNumGroups(), 0.0); + vector temp; + for (int j = 0; j < m->getNumGroups(); j++) { temp.push_back(0.0); } + //stddists[i].resize(m->getNumGroups(), 0.0); + stddists.push_back(temp); } + //flip it so you can print it int count = 0; @@ -509,13 +500,18 @@ int UnifracWeightedCommand::getAverageSTDMatrices(vector< vector >& dist } } - string aveFileName = outputDir + m->getSimpleName(treefile) + toString(treeNum+1) + ".weighted.ave." + getOutputFileNameTag("phylip"); + map variables; + variables["[filename]"] = outputDir + m->getSimpleName(treefile); + variables["[tag]"] = toString(treeNum+1); + variables["[tag2]"] = "weighted.ave"; + string aveFileName = getOutputFileName("phylip",variables); if (outputForm != "column") { outputNames.push_back(aveFileName); outputTypes["phylip"].push_back(aveFileName); } else { outputNames.push_back(aveFileName); outputTypes["column"].push_back(aveFileName); } ofstream out; m->openOutputFile(aveFileName, out); - string stdFileName = outputDir + m->getSimpleName(treefile) + toString(treeNum+1) + ".weighted.std." + getOutputFileNameTag("phylip"); + variables["[tag2]"] = "weighted.std"; + string stdFileName = getOutputFileName("phylip",variables); if (outputForm != "column") { outputNames.push_back(stdFileName); outputTypes["phylip"].push_back(stdFileName); } else { outputNames.push_back(stdFileName); outputTypes["column"].push_back(stdFileName); } ofstream outStd; @@ -606,7 +602,11 @@ int UnifracWeightedCommand::getConsensusTrees(vector< vector >& dists, i Tree* conTree = con.getTree(newTrees); //create a new filename - string conFile = outputDir + m->getRootName(m->getSimpleName(treefile)) + toString(treeNum+1) + ".weighted.cons." + getOutputFileNameTag("tree"); + map variables; + variables["[filename]"] = outputDir + m->getRootName(m->getSimpleName(treefile)); + variables["[tag]"] = toString(treeNum+1); + variables["[tag2]"] = "weighted.cons"; + string conFile = getOutputFileName("tree",variables); outputNames.push_back(conFile); outputTypes["tree"].push_back(conFile); ofstream outTree; m->openOutputFile(conFile, outTree); @@ -630,7 +630,11 @@ vector UnifracWeightedCommand::buildTrees(vector< vector >& dists vector trees; //create a new filename - string outputFile = outputDir + m->getRootName(m->getSimpleName(treefile)) + toString(treeNum+1) + ".weighted.all." + getOutputFileNameTag("tree"); + map variables; + variables["[filename]"] = outputDir + m->getRootName(m->getSimpleName(treefile)); + variables["[tag]"] = toString(treeNum+1); + variables["[tag2]"] = "weighted.all"; + string outputFile = getOutputFileName("tree",variables); outputNames.push_back(outputFile); outputTypes["tree"].push_back(outputFile); ofstream outAll; @@ -694,39 +698,32 @@ int UnifracWeightedCommand::runRandomCalcs(Tree* thisTree, vector usersS lines.clear(); -#if defined (__APPLE__) || (__MACH__) || (linux) || (__linux) || (__linux__) || (__unix__) || (__unix) - if(processors != 1){ - int numPairs = namesOfGroupCombos.size(); - int numPairsPerProcessor = numPairs / processors; + //breakdown work between processors + 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)); - } + 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)); } -#endif //get scores for random trees for (int j = 0; j < iters; j++) { - cout << j << endl; -#if defined (__APPLE__) || (__MACH__) || (linux) || (__linux) || (__linux__) || (__unix__) || (__unix) - if(processors == 1){ - driver(thisTree, namesOfGroupCombos, 0, namesOfGroupCombos.size(), rScores); - }else{ +//#if defined (__APPLE__) || (__MACH__) || (linux) || (__linux) || (__linux__) || (__unix__) || (__unix) + //if(processors == 1){ + // driver(thisTree, namesOfGroupCombos, 0, namesOfGroupCombos.size(), rScores); + // }else{ createProcesses(thisTree, namesOfGroupCombos, rScores); - } -#else - driver(thisTree, namesOfGroupCombos, 0, namesOfGroupCombos.size(), rScores); -#endif + // } +//#else + //driver(thisTree, namesOfGroupCombos, 0, namesOfGroupCombos.size(), rScores); +//#endif + if (m->control_pressed) { delete ct; for (int i = 0; i < T.size(); i++) { delete T[i]; } delete output; outSum.close(); for (int i = 0; i < outputNames.size(); i++) { m->mothurRemove(outputNames[i]); } return 0; } - //report progress - // m->mothurOut("Iter: " + toString(j+1)); m->mothurOutEndLine(); } lines.clear(); @@ -762,12 +759,11 @@ int UnifracWeightedCommand::runRandomCalcs(Tree* thisTree, vector usersS int UnifracWeightedCommand::createProcesses(Tree* t, vector< vector > namesOfGroupCombos, vector< vector >& scores) { try { -#if defined (__APPLE__) || (__MACH__) || (linux) || (__linux) || (__linux__) || (__unix__) || (__unix) - int process = 1; + int process = 1; vector 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(); @@ -813,9 +809,53 @@ int UnifracWeightedCommand::createProcesses(Tree* t, vector< vector > na in.close(); m->mothurRemove(s); } +#else + //fill in functions + vector pDataArray; + DWORD dwThreadIdArray[processors-1]; + HANDLE hThreadArray[processors-1]; + vector cts; + vector trees; - return 0; -#endif + //Create processor worker threads. + for( int i=1; icopy(ct); + Tree* copyTree = new Tree(copyCount); + copyTree->getCopy(t); + + cts.push_back(copyCount); + trees.push_back(copyTree); + + vector< vector > copyScores = rScores; + + weightedRandomData* tempweighted = new weightedRandomData(m, lines[i].start, lines[i].num, namesOfGroupCombos, copyTree, copyCount, includeRoot, copyScores); + pDataArray.push_back(tempweighted); + processIDS.push_back(i); + + hThreadArray[i-1] = CreateThread(NULL, 0, MyWeightedRandomThreadFunction, pDataArray[i-1], 0, &dwThreadIdArray[i-1]); + } + + driver(t, namesOfGroupCombos, lines[0].start, lines[0].num, scores); + + //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 = pDataArray[i]->start; j < (pDataArray[i]->start+pDataArray[i]->num); j++) { + scores[j].push_back(pDataArray[i]->scores[j][pDataArray[i]->scores[j].size()-1]); + } + delete cts[i]; + delete trees[i]; + CloseHandle(hThreadArray[i]); + delete pDataArray[i]; + } + + +#endif + + return 0; } catch(exception& e) { m->errorOut(e, "UnifracWeightedCommand", "createProcesses"); @@ -938,14 +978,20 @@ void UnifracWeightedCommand::createPhylipFile() { //for each tree for (int i = 0; i < T.size(); i++) { - string phylipFileName; - if ((outputForm == "lt") || (outputForm == "square")) { - phylipFileName = outputDir + m->getSimpleName(treefile) + toString(i+1) + ".weighted.phylip." + getOutputFileNameTag("phylip"); - outputNames.push_back(phylipFileName); outputTypes["phylip"].push_back(phylipFileName); - }else { //column - phylipFileName = outputDir + m->getSimpleName(treefile) + toString(i+1) + ".weighted.column." + getOutputFileNameTag("column"); - outputNames.push_back(phylipFileName); outputTypes["column"].push_back(phylipFileName); - } + string phylipFileName; + map variables; + variables["[filename]"] = outputDir + m->getSimpleName(treefile); + variables["[tag]"] = toString(i+1); + if ((outputForm == "lt") || (outputForm == "square")) { + variables["[tag2]"] = "weighted.phylip"; + phylipFileName = getOutputFileName("phylip",variables); + outputNames.push_back(phylipFileName); outputTypes["phylip"].push_back(phylipFileName); + }else { //column + variables["[tag2]"] = "weighted.column"; + phylipFileName = getOutputFileName("column",variables); + outputNames.push_back(phylipFileName); outputTypes["column"].push_back(phylipFileName); + } + ofstream out; m->openOutputFile(phylipFileName, out);