2 // getmetacommunitycommand.cpp
5 // Created by SarahsWork on 4/9/13.
6 // Copyright (c) 2013 Schloss Lab. All rights reserved.
9 #include "getmetacommunitycommand.h"
12 //**********************************************************************************************************************
13 vector<string> GetMetaCommunityCommand::setParameters(){
15 CommandParameter pshared("shared", "InputTypes", "", "", "none", "none", "none","outputType",false,true); parameters.push_back(pshared);
16 CommandParameter pgroups("groups", "String", "", "", "", "", "","",false,false); parameters.push_back(pgroups);
17 CommandParameter plabel("label", "String", "", "", "", "", "","",false,false); parameters.push_back(plabel);
18 CommandParameter pminpartitions("minpartitions", "Number", "", "5", "", "", "","",false,false,true); parameters.push_back(pminpartitions);
19 CommandParameter pmaxpartitions("maxpartitions", "Number", "", "100", "", "", "","",false,false,true); parameters.push_back(pmaxpartitions);
20 CommandParameter poptimizegap("optimizegap", "Number", "", "3", "", "", "","",false,false,true); parameters.push_back(poptimizegap);
21 CommandParameter pprocessors("processors", "Number", "", "1", "", "", "","",false,false,true); parameters.push_back(pprocessors);
22 CommandParameter pinputdir("inputdir", "String", "", "", "", "", "","",false,false); parameters.push_back(pinputdir);
23 CommandParameter poutputdir("outputdir", "String", "", "", "", "", "","",false,false); parameters.push_back(poutputdir);
25 vector<string> myArray;
26 for (int i = 0; i < parameters.size(); i++) { myArray.push_back(parameters[i].name); }
30 m->errorOut(e, "NewCommand", "setParameters");
34 //**********************************************************************************************************************
35 string GetMetaCommunityCommand::getHelpString(){
37 string helpString = "";
38 helpString += "The get.communitytype command parameters are shared, label, groups, minpartitions, maxpartitions, optimizegap and processors. The shared file is required. \n";
39 helpString += "The label parameter is used to analyze specific labels in your input. labels are separated by dashes.\n";
40 helpString += "The groups parameter allows you to specify which of the groups in your shared file you would like analyzed. Group names are separated by dashes.\n";
41 helpString += "The minpartitions parameter is used to .... Default=5.\n";
42 helpString += "The maxpartitions parameter is used to .... Default=10.\n";
43 helpString += "The optimizegap parameter is used to .... Default=3.\n";
44 helpString += "The processors parameter allows you to specify number of processors to use. The default is 1.\n";
45 helpString += "The get.communitytype command should be in the following format: get.communitytype(shared=yourSharedFile).\n";
49 m->errorOut(e, "GetMetaCommunityCommand", "getHelpString");
53 //**********************************************************************************************************************
54 string GetMetaCommunityCommand::getOutputPattern(string type) {
58 if (type == "fit") { pattern = "[filename],[distance],mix.fit"; }
59 else if (type == "relabund") { pattern = "[filename],[distance],[tag],mix.relabund"; }
60 else if (type == "design") { pattern = "[filename],[distance],mix.design"; }
61 else if (type == "matrix") { pattern = "[filename],[distance],[tag],mix.posterior"; }
62 else if (type == "parameters") { pattern = "[filename],[distance],mix.parameters"; }
63 else if (type == "summary") { pattern = "[filename],[distance],mix.summary"; }
64 else { m->mothurOut("[ERROR]: No definition for type " + type + " output pattern.\n"); m->control_pressed = true; }
69 m->errorOut(e, "GetMetaCommunityCommand", "getOutputPattern");
73 //**********************************************************************************************************************
74 GetMetaCommunityCommand::GetMetaCommunityCommand(){
76 abort = true; calledHelp = true;
78 vector<string> tempOutNames;
79 outputTypes["fit"] = tempOutNames;
80 outputTypes["relabund"] = tempOutNames;
81 outputTypes["matrix"] = tempOutNames;
82 outputTypes["design"] = tempOutNames;
83 outputTypes["parameters"] = tempOutNames;
84 outputTypes["summary"] = tempOutNames;
87 m->errorOut(e, "GetMetaCommunityCommand", "GetMetaCommunityCommand");
91 //**********************************************************************************************************************
92 GetMetaCommunityCommand::GetMetaCommunityCommand(string option) {
94 abort = false; calledHelp = false;
97 //allow user to run help
98 if(option == "help") { help(); abort = true; calledHelp = true; }
99 else if(option == "citation") { citation(); abort = true; calledHelp = true;}
102 //valid paramters for this command
103 vector<string> myArray = setParameters();
105 OptionParser parser(option);
106 map<string,string> parameters = parser.getParameters();
108 ValidParameters validParameter;
109 map<string,string>::iterator it;
110 //check to make sure all parameters are valid for command
111 for (it = parameters.begin(); it != parameters.end(); it++) {
112 if (validParameter.isValidParameter(it->first, myArray, it->second) != true) { abort = true; }
115 vector<string> tempOutNames;
116 outputTypes["fit"] = tempOutNames;
117 outputTypes["relabund"] = tempOutNames;
118 outputTypes["matrix"] = tempOutNames;
119 outputTypes["design"] = tempOutNames;
120 outputTypes["parameters"] = tempOutNames;
121 outputTypes["summary"] = tempOutNames;
123 //if the user changes the input directory command factory will send this info to us in the output parameter
124 string inputDir = validParameter.validFile(parameters, "inputdir", false);
125 if (inputDir == "not found"){ inputDir = ""; }
128 it = parameters.find("shared");
129 if(it != parameters.end()){
130 path = m->hasPath(it->second);
131 if (path == "") { parameters["shared"] = inputDir + it->second; }
135 //get shared file, it is required
136 sharedfile = validParameter.validFile(parameters, "shared", true);
137 if (sharedfile == "not open") { sharedfile = ""; abort = true; }
138 else if (sharedfile == "not found") {
139 //if there is a current shared file, use it
140 sharedfile = m->getSharedFile();
141 if (sharedfile != "") { m->mothurOut("Using " + sharedfile + " as input file for the shared parameter."); m->mothurOutEndLine(); }
142 else { m->mothurOut("You have no current sharedfile and the shared parameter is required."); m->mothurOutEndLine(); abort = true; }
143 }else { m->setSharedFile(sharedfile); }
145 //if the user changes the output directory command factory will send this info to us in the output parameter
146 outputDir = validParameter.validFile(parameters, "outputdir", false); if (outputDir == "not found"){
147 outputDir = m->hasPath(sharedfile); //if user entered a file with a path then preserve it
150 string temp = validParameter.validFile(parameters, "minpartitions", false); if (temp == "not found"){ temp = "5"; }
151 m->mothurConvert(temp, minpartitions);
153 temp = validParameter.validFile(parameters, "maxpartitions", false); if (temp == "not found"){ temp = "10"; }
154 m->mothurConvert(temp, maxpartitions);
156 temp = validParameter.validFile(parameters, "optimizegap", false); if (temp == "not found"){ temp = "3"; }
157 m->mothurConvert(temp, optimizegap);
159 temp = validParameter.validFile(parameters, "processors", false); if (temp == "not found"){ temp = m->getProcessors(); }
160 m->setProcessors(temp);
161 m->mothurConvert(temp, processors);
163 string groups = validParameter.validFile(parameters, "groups", false);
164 if (groups == "not found") { groups = ""; }
165 else { m->splitAtDash(groups, Groups); }
166 m->setGroups(Groups);
168 string label = validParameter.validFile(parameters, "label", false);
169 if (label == "not found") { label = ""; }
171 if(label != "all") { m->splitAtDash(label, labels); allLines = 0; }
172 else { allLines = 1; }
177 catch(exception& e) {
178 m->errorOut(e, "GetMetaCommunityCommand", "GetMetaCommunityCommand");
182 //**********************************************************************************************************************
184 int GetMetaCommunityCommand::execute(){
187 if (abort == true) { if (calledHelp) { return 0; } return 2; }
189 InputData input(sharedfile, "sharedfile");
190 vector<SharedRAbundVector*> lookup = input.getSharedRAbundVectors();
191 string lastLabel = lookup[0]->getLabel();
193 //if the users enters label "0.06" and there is no "0.06" in their file use the next lowest label.
194 set<string> processedLabels;
195 set<string> userLabels = labels;
197 //as long as you are not at the end of the file or done wih the lines you want
198 while((lookup[0] != NULL) && ((allLines == 1) || (userLabels.size() != 0))) {
200 if (m->control_pressed) { for (int i = 0; i < lookup.size(); i++) { delete lookup[i]; } return 0; }
202 if(allLines == 1 || labels.count(lookup[0]->getLabel()) == 1){
204 m->mothurOut(lookup[0]->getLabel()); m->mothurOutEndLine();
206 createProcesses(lookup);
208 processedLabels.insert(lookup[0]->getLabel());
209 userLabels.erase(lookup[0]->getLabel());
212 if ((m->anyLabelsToProcess(lookup[0]->getLabel(), userLabels, "") == true) && (processedLabels.count(lastLabel) != 1)) {
213 string saveLabel = lookup[0]->getLabel();
215 for (int i = 0; i < lookup.size(); i++) { delete lookup[i]; }
216 lookup = input.getSharedRAbundVectors(lastLabel);
217 m->mothurOut(lookup[0]->getLabel()); m->mothurOutEndLine();
219 createProcesses(lookup);
221 processedLabels.insert(lookup[0]->getLabel());
222 userLabels.erase(lookup[0]->getLabel());
224 //restore real lastlabel to save below
225 lookup[0]->setLabel(saveLabel);
228 lastLabel = lookup[0]->getLabel();
229 //prevent memory leak
230 for (int i = 0; i < lookup.size(); i++) { delete lookup[i]; lookup[i] = NULL; }
232 if (m->control_pressed) { return 0; }
234 //get next line to process
235 lookup = input.getSharedRAbundVectors();
238 if (m->control_pressed) { return 0; }
240 //output error messages about any remaining user labels
241 set<string>::iterator it;
242 bool needToRun = false;
243 for (it = userLabels.begin(); it != userLabels.end(); it++) {
244 m->mothurOut("Your file does not include the label " + *it);
245 if (processedLabels.count(lastLabel) != 1) {
246 m->mothurOut(". I will use " + lastLabel + "."); m->mothurOutEndLine();
249 m->mothurOut(". Please refer to " + lastLabel + "."); m->mothurOutEndLine();
253 //run last label if you need to
254 if (needToRun == true) {
255 for (int i = 0; i < lookup.size(); i++) { if (lookup[i] != NULL) { delete lookup[i]; } }
256 lookup = input.getSharedRAbundVectors(lastLabel);
258 m->mothurOut(lookup[0]->getLabel()); m->mothurOutEndLine();
260 createProcesses(lookup);
262 for (int i = 0; i < lookup.size(); i++) { delete lookup[i]; }
265 //output files created by command
266 m->mothurOutEndLine();
267 m->mothurOut("Output File Names: "); m->mothurOutEndLine();
268 for (int i = 0; i < outputNames.size(); i++) { m->mothurOut(outputNames[i]); m->mothurOutEndLine(); }
269 m->mothurOutEndLine();
273 catch(exception& e) {
274 m->errorOut(e, "GetMetaCommunityCommand", "execute");
278 //**********************************************************************************************************************
279 int GetMetaCommunityCommand::createProcesses(vector<SharedRAbundVector*>& thislookup){
282 #if defined (__APPLE__) || (__MACH__) || (linux) || (__linux) || (__linux__) || (__unix__) || (__unix)
284 processors=1; //qFinderDMM not thread safe
287 vector<int> processIDS;
290 int minPartition = 0;
293 if (maxpartitions < processors) { processors = maxpartitions; }
295 map<string, string> variables;
296 variables["[filename]"] = outputDir + m->getRootName(m->getSimpleName(sharedfile));
297 variables["[distance]"] = thislookup[0]->getLabel();
298 string outputFileName = getOutputFileName("fit", variables);
299 outputNames.push_back(outputFileName); outputTypes["fit"].push_back(outputFileName);
301 //divide the partitions between the processors
302 vector< vector<int> > dividedPartitions;
303 vector< vector<string> > rels, matrix;
304 vector<string> doneFlags;
305 dividedPartitions.resize(processors);
306 rels.resize(processors);
307 matrix.resize(processors);
309 //for each file group figure out which process will complete it
310 //want to divide the load intelligently so the big files are spread between processes
311 for (int i=1; i<=maxpartitions; i++) {
313 int processToAssign = (i+1) % processors;
314 if (processToAssign == 0) { processToAssign = processors; }
316 if (m->debug) { m->mothurOut("[DEBUG]: assigning " + toString(i) + " to process " + toString(processToAssign-1) + "\n"); }
317 dividedPartitions[(processToAssign-1)].push_back(i);
319 variables["[tag]"] = toString(i);
320 string relName = getOutputFileName("relabund", variables);
321 string mName = getOutputFileName("matrix", variables);
322 rels[(processToAssign-1)].push_back(relName);
323 matrix[(processToAssign-1)].push_back(mName);
326 for (int i = 0; i < processors; i++) { //read from everyone elses, just write to yours
327 string tempDoneFile = m->getRootName(m->getSimpleName(sharedfile)) + toString(i) + ".done.temp";
328 doneFlags.push_back(tempDoneFile);
330 m->openOutputFile(tempDoneFile, out); //clear out
335 #if defined (__APPLE__) || (__MACH__) || (linux) || (__linux) || (__linux__) || (__unix__) || (__unix)
337 //loop through and create all the processes you want
338 while (process != processors) {
342 processIDS.push_back(pid); //create map from line number to pid so you can append files in correct order later
346 num = processDriver(thislookup, dividedPartitions[process], (outputFileName + toString(getpid())), rels[process], matrix[process], doneFlags, process);
348 //pass numSeqs to parent
350 string tempFile = toString(getpid()) + ".outputNames.temp";
351 m->openOutputFile(tempFile, out);
353 out << outputNames.size() << endl;
354 for (int i = 0; i < outputNames.size(); i++) { out << outputNames[i] << endl; }
359 m->mothurOut("[ERROR]: unable to spawn the necessary processes."); m->mothurOutEndLine();
360 for (int i = 0; i < processIDS.size(); i++) { kill (processIDS[i], SIGINT); }
366 m->mothurOut("K\tNLE\t\tlogDet\tBIC\t\tAIC\t\tLaplace\n");
367 minPartition = processDriver(thislookup, dividedPartitions[0], outputFileName, rels[0], matrix[0], doneFlags, 0);
369 //force parent to wait until all the processes are done
370 for (int i=0;i<processIDS.size();i++) {
371 int temp = processIDS[i];
375 vector<string> tempOutputNames = outputNames;
376 for (int i=0;i<processIDS.size();i++) {
378 string tempFile = toString(processIDS[i]) + ".outputNames.temp";
379 m->openInputFile(tempFile, in);
382 in >> tempNum; m->gobble(in);
383 if (tempNum < minPartition) { minPartition = tempNum; }
384 in >> tempNum; m->gobble(in);
385 for (int i = 0; i < tempNum; i++) {
386 string tempName = "";
387 in >> tempName; m->gobble(in);
388 tempOutputNames.push_back(tempName);
391 in.close(); m->mothurRemove(tempFile);
393 m->appendFilesWithoutHeaders(outputFileName + toString(processIDS[i]), outputFileName);
394 m->mothurRemove(outputFileName + toString(processIDS[i]));
397 if (processors > 1) {
399 for (int i = 0; i < tempOutputNames.size(); i++) { //remove files if needed
400 string name = tempOutputNames[i];
401 vector<string> parts;
402 m->splitAtChar(name, parts, '.');
404 if (((parts[parts.size()-1] == "relabund") || (parts[parts.size()-1] == "posterior")) && (parts[parts.size()-2] == "mix")) {
405 string tempNum = parts[parts.size()-3];
406 int num; m->mothurConvert(tempNum, num);
407 //if (num > minPartition) {
408 // m->mothurRemove(tempOutputNames[i]);
409 // keep = false; if (m->debug) { m->mothurOut("[DEBUG]: removing " + tempOutputNames[i] + ".\n"); }
412 if (keep) { outputNames.push_back(tempOutputNames[i]); }
417 m->openInputFile(outputFileName, in);
418 string headers = m->getline(in); m->gobble(in);
420 map<int, string> file;
422 string numString, line;
424 in >> numString; line = m->getline(in); m->gobble(in);
425 m->mothurConvert(numString, num);
430 m->openOutputFile(outputFileName, out);
431 out << headers << endl;
432 for (map<int, string>::iterator it = file.begin(); it != file.end(); it++) {
433 out << it->first << '\t' << it->second << endl;
434 if (m->debug) { m->mothurOut("[DEBUG]: printing: " + toString(it->first) + '\t' + it->second + ".\n"); }
441 vector<metaCommunityData*> pDataArray;
442 DWORD dwThreadIdArray[processors-1];
443 HANDLE hThreadArray[processors-1];
445 //Create processor worker threads.
446 for( int i=1; i<processors; i++ ){
448 //make copy of lookup so we don't get access violations
449 vector<SharedRAbundVector*> newLookup;
450 for (int k = 0; k < thislookup.size(); k++) {
451 SharedRAbundVector* temp = new SharedRAbundVector();
452 temp->setLabel(thislookup[k]->getLabel());
453 temp->setGroup(thislookup[k]->getGroup());
454 newLookup.push_back(temp);
458 for (int k = 0; k < thislookup[0]->getNumBins(); k++) {
459 if (m->control_pressed) { for (int j = 0; j < newLookup.size(); j++) { delete newLookup[j]; } return 0; }
460 for (int j = 0; j < thislookup.size(); j++) { newLookup[j]->push_back(thislookup[j]->getAbundance(k), thislookup[j]->getGroup()); }
463 processIDS.push_back(i);
465 // Allocate memory for thread data.
466 metaCommunityData* tempMeta = new metaCommunityData(newLookup, m, dividedPartitions[i], outputFileName + toString(i), rels[i], matrix[i], minpartitions, maxpartitions, optimizegap);
467 pDataArray.push_back(tempMeta);
469 hThreadArray[i] = CreateThread(NULL, 0, MyMetaCommunityThreadFunction, pDataArray[i], 0, &dwThreadIdArray[i]);
473 minPartition = processDriver(thislookup, dividedPartitions[0], outputFileName, rels[0], matrix[0]);
475 //Wait until all threads have terminated.
476 WaitForMultipleObjects(processors-1, hThreadArray, TRUE, INFINITE);
478 //Close all thread handles and free memory allocations.
479 for(int i=0; i < pDataArray.size(); i++){
480 if (pDataArray[i]->minPartition < minPartition) { minPartition = pDataArray[i]->minPartition; }
481 for (int j = 0; j < pDataArray[i]->outputNames.size(); j++) {
482 outputNames.push_back(pDataArray[i]->outputNames[j]);
484 m->appendFilesWithoutHeaders(outputFileName + toString(processIDS[i]), outputFileName);
485 m->mothurRemove(outputFileName + toString(processIDS[i]));
486 CloseHandle(hThreadArray[i]);
487 delete pDataArray[i];
490 m->mothurOut("K\tNLE\t\tlogDet\tBIC\t\tAIC\t\tLaplace\n");
491 minPartition = processDriver(thislookup, dividedPartitions[0], outputFileName, rels[0], matrix[0], doneFlags, 0);
493 for (int i = 0; i < processors; i++) { //read from everyone elses, just write to yours
494 string tempDoneFile = m->getRootName(m->getSimpleName(sharedfile)) + toString(i) + ".done.temp";
495 m->mothurRemove(tempDoneFile);
498 if (m->control_pressed) { return 0; }
500 if (m->debug) { m->mothurOut("[DEBUG]: minPartition = " + toString(minPartition) + "\n"); }
502 //run generate Summary function for smallest minPartition
503 variables["[tag]"] = toString(minPartition);
504 generateSummaryFile(minPartition, variables);
509 catch(exception& e) {
510 m->errorOut(e, "GetMetaCommunityCommand", "createProcesses");
514 //**********************************************************************************************************************
515 int GetMetaCommunityCommand::processDriver(vector<SharedRAbundVector*>& thislookup, vector<int>& parts, string outputFileName, vector<string> relabunds, vector<string> matrix, vector<string> doneFlags, int processID){
518 double minLaplace = 1e10;
519 int minPartition = 0;
522 m->openOutputFile(outputFileName, fitData);
523 fitData.setf(ios::fixed, ios::floatfield);
524 fitData.setf(ios::showpoint);
525 cout.setf(ios::fixed, ios::floatfield);
526 cout.setf(ios::showpoint);
528 vector< vector<int> > sharedMatrix;
529 vector<string> thisGroups;
530 for (int i = 0; i < thislookup.size(); i++) { sharedMatrix.push_back(thislookup[i]->getAbundances()); thisGroups.push_back(thislookup[i]->getGroup()); }
532 fitData << "K\tNLE\tlogDet\tBIC\tAIC\tLaplace" << endl;
534 for(int i=0;i<parts.size();i++){
536 int numPartitions = parts[i];
538 if (m->debug) { m->mothurOut("[DEBUG]: running partition " + toString(numPartitions) + "\n"); }
540 if (m->control_pressed) { break; }
542 //check to see if anyone else is done
543 for (int j = 0; j < doneFlags.size(); j++) {
544 if (!m->isBlank(doneFlags[j])) { //another process has finished
545 //are they done at a lower partition?
547 m->openInputFile(doneFlags[j], in);
549 in >> tempNum; in.close();
550 if (tempNum < numPartitions) { break; } //quit, because someone else has finished
554 qFinderDMM findQ(sharedMatrix, numPartitions);
556 double laplace = findQ.getLaplace();
557 m->mothurOut(toString(numPartitions) + '\t');
558 cout << setprecision (2) << findQ.getNLL() << '\t' << findQ.getLogDet() << '\t';
559 m->mothurOutJustToLog(toString(findQ.getNLL()) + '\t' + toString(findQ.getLogDet()) + '\t');
560 cout << findQ.getBIC() << '\t' << findQ.getAIC() << '\t' << laplace;
561 m->mothurOutJustToLog(toString(findQ.getBIC()) + '\t' + toString(findQ.getAIC()) + '\t' + toString(laplace));
563 fitData << numPartitions << '\t';
564 fitData << setprecision (2) << findQ.getNLL() << '\t' << findQ.getLogDet() << '\t';
565 fitData << findQ.getBIC() << '\t' << findQ.getAIC() << '\t' << laplace << endl;
567 if(laplace < minLaplace){
568 minPartition = numPartitions;
569 minLaplace = laplace;
572 m->mothurOutEndLine();
574 string relabund = relabunds[i];
575 outputNames.push_back(relabund); outputTypes["relabund"].push_back(relabund);
576 string matrixName = matrix[i];
577 outputNames.push_back(matrixName); outputTypes["matrix"].push_back(matrixName);
579 findQ.printZMatrix(matrixName, thisGroups);
580 findQ.printRelAbund(relabund, m->currentBinLabels);
582 if(optimizegap != -1 && (numPartitions - minPartition) >= optimizegap && numPartitions >= minpartitions){
583 string tempDoneFile = m->getRootName(m->getSimpleName(sharedfile)) + toString(processID) + ".done.temp";
585 m->openOutputFile(tempDoneFile, outDone);
586 outDone << minPartition << endl;
595 if (m->control_pressed) { return 0; }
599 catch(exception& e) {
600 m->errorOut(e, "GetMetaCommunityCommand", "processDriver");
604 /**************************************************************************************************/
606 vector<double> GetMetaCommunityCommand::generateDesignFile(int numPartitions, map<string,string> variables){
608 vector<double> piValues(numPartitions, 0);
611 variables["[tag]"] = toString(numPartitions);
612 string input = getOutputFileName("matrix", variables);
613 m->openInputFile(input, postFile);//((fileRoot + toString(numPartitions) + "mix.posterior").c_str()); //matrix file
614 variables.erase("[tag]");
615 string outputFileName = getOutputFileName("design", variables);
617 m->openOutputFile(outputFileName, designFile);
618 outputNames.push_back(outputFileName); outputTypes["design"].push_back(outputFileName);
621 vector<string> titles(numPartitions);
623 for(int i=0;i<numPartitions;i++){ postFile >> titles[i]; }
631 if (m->control_pressed) { break; }
633 double maxPosterior = 0.0000;
634 int maxPartition = -1;
636 postFile >> sampleName;
638 for(int i=0;i<numPartitions;i++){
640 postFile >> posterior;
641 if(posterior > maxPosterior){
642 maxPosterior = posterior;
645 piValues[i] += posterior;
649 designFile << sampleName << '\t' << titles[maxPartition] << endl;
654 for(int i=0;i<numPartitions;i++){
655 piValues[i] /= (double)numSamples;
664 catch(exception& e) {
665 m->errorOut(e, "GetMetaCommunityCommand", "generateDesignFile");
670 /**************************************************************************************************/
672 inline bool summaryFunction(summaryData i, summaryData j){ return i.difference > j.difference; }
674 /**************************************************************************************************/
675 int GetMetaCommunityCommand::generateSummaryFile(int numPartitions, map<string,string> v){
677 vector<summaryData> summary;
679 vector<double> pMean(numPartitions, 0);
680 vector<double> pLCI(numPartitions, 0);
681 vector<double> pUCI(numPartitions, 0);
684 double mean, lci, uci;
687 vector<double> piValues = generateDesignFile(numPartitions, v);
689 ifstream referenceFile;
690 map<string, string> variables;
691 variables["[filename]"] = v["[filename]"];
692 variables["[distance]"] = v["[distance]"];
693 variables["[tag]"] = "1";
694 string reference = getOutputFileName("relabund", variables);
695 m->openInputFile(reference, referenceFile); //((fileRoot + label + ".1mix.relabund").c_str());
696 variables["[tag]"] = toString(numPartitions);
697 string partFile = getOutputFileName("relabund", variables);
698 ifstream partitionFile;
699 m->openInputFile(partFile, partitionFile); //((fileRoot + toString(numPartitions) + "mix.relabund").c_str());
701 header = m->getline(referenceFile);
702 header = m->getline(partitionFile);
703 stringstream head(header);
706 vector<string> thetaValues(numPartitions, "");
707 for(int i=0;i<numPartitions;i++){
708 head >> label >> dummy >> dummy;
709 thetaValues[i] = label.substr(label.find_last_of('_')+1);
713 vector<double> partitionDiff(numPartitions, 0.0000);
715 while(referenceFile){
717 if (m->control_pressed) { break; }
719 referenceFile >> name >> mean >> lci >> uci;
721 summaryData tempData;
722 tempData.name = name;
723 tempData.refMean = mean;
725 double difference = 0.0000;
727 partitionFile >> name;
728 for(int j=0;j<numPartitions;j++){
729 partitionFile >> pMean[j] >> pLCI[j] >> pUCI[j];
730 difference += abs(mean - pMean[j]);
731 partitionDiff[j] += abs(mean - pMean[j]);;
734 tempData.partMean = pMean;
735 tempData.partLCI = pLCI;
736 tempData.partUCI = pUCI;
737 tempData.difference = difference;
738 summary.push_back(tempData);
740 m->gobble(referenceFile);
741 m->gobble(partitionFile);
743 referenceFile.close();
744 partitionFile.close();
746 if (m->control_pressed) { return 0; }
748 int numOTUs = (int)summary.size();
750 sort(summary.begin(), summary.end(), summaryFunction);
752 variables.erase("[tag]");
753 string outputFileName = getOutputFileName("parameters", variables);
754 outputNames.push_back(outputFileName); outputTypes["parameters"].push_back(outputFileName);
756 ofstream parameterFile;
757 m->openOutputFile(outputFileName, parameterFile); //((fileRoot + "mix.parameters").c_str());
758 parameterFile.setf(ios::fixed, ios::floatfield);
759 parameterFile.setf(ios::showpoint);
761 double totalDifference = 0.0000;
762 parameterFile << "Part\tDif2Ref_i\ttheta_i\tpi_i\n";
763 for(int i=0;i<numPartitions;i++){
764 if (m->control_pressed) { break; }
765 parameterFile << i+1 << '\t' << setprecision(2) << partitionDiff[i] << '\t' << thetaValues[i] << '\t' << piValues[i] << endl;
766 totalDifference += partitionDiff[i];
768 parameterFile.close();
770 if (m->control_pressed) { return 0; }
772 string summaryFileName = getOutputFileName("summary", variables);
773 outputNames.push_back(summaryFileName); outputTypes["summary"].push_back(summaryFileName);
775 ofstream summaryFile;
776 m->openOutputFile(summaryFileName, summaryFile); //((fileRoot + "mix.summary").c_str());
777 summaryFile.setf(ios::fixed, ios::floatfield);
778 summaryFile.setf(ios::showpoint);
781 summaryFile << "OTU\tP0.mean";
782 for(int i=0;i<numPartitions;i++){
783 summaryFile << "\tP" << i+1 << ".mean\tP" << i+1 << ".lci\tP" << i+1 << ".uci";
785 summaryFile << "\tDifference\tCumFraction" << endl;
787 double cumDiff = 0.0000;
789 for(int i=0;i<numOTUs;i++){
790 if (m->control_pressed) { break; }
791 summaryFile << summary[i].name << setprecision(2) << '\t' << summary[i].refMean;
792 for(int j=0;j<numPartitions;j++){
793 summaryFile << '\t' << summary[i].partMean[j] << '\t' << summary[i].partLCI[j] << '\t' << summary[i].partUCI[j];
796 cumDiff += summary[i].difference/totalDifference;
797 summaryFile << '\t' << summary[i].difference << '\t' << cumDiff << endl;
804 catch(exception& e) {
805 m->errorOut(e, "GetMetaCommunityCommand", "generateSummaryFile");
810 //**********************************************************************************************************************