getmetacommunitycommand.cpp

   1 //
   2 //  getmetacommunitycommand.cpp
   3 //  Mothur
   4 //
   5 //  Created by SarahsWork on 4/9/13.
   6 //  Copyright (c) 2013 Schloss Lab. All rights reserved.
   7 //
   8
   9 #include "getmetacommunitycommand.h"
  10
  11
  12 //**********************************************************************************************************************
  13 vector<string> GetMetaCommunityCommand::setParameters(){
  14         try {
  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);
  24
  25                 vector<string> myArray;
  26                 for (int i = 0; i < parameters.size(); i++) {   myArray.push_back(parameters[i].name);          }
  27                 return myArray;
  28         }
  29         catch(exception& e) {
  30                 m->errorOut(e, "NewCommand", "setParameters");
  31                 exit(1);
  32         }
  33 }
  34 //**********************************************************************************************************************
  35 string GetMetaCommunityCommand::getHelpString(){
  36         try {
  37                 string helpString = "";
  38                 helpString += "The get.metacommunity 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.metacommunity command should be in the following format: get.metacommunity(shared=yourSharedFile).\n";
  46                 return helpString;
  47         }
  48         catch(exception& e) {
  49                 m->errorOut(e, "GetMetaCommunityCommand", "getHelpString");
  50                 exit(1);
  51         }
  52 }
  53 //**********************************************************************************************************************
  54 string GetMetaCommunityCommand::getOutputPattern(string type) {
  55     try {
  56         string pattern = "";
  57
  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;  }
  65
  66         return pattern;
  67     }
  68     catch(exception& e) {
  69         m->errorOut(e, "GetMetaCommunityCommand", "getOutputPattern");
  70         exit(1);
  71     }
  72 }
  73 //**********************************************************************************************************************
  74 GetMetaCommunityCommand::GetMetaCommunityCommand(){
  75         try {
  76                 abort = true; calledHelp = true;
  77                 setParameters();
  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;
  85         }
  86         catch(exception& e) {
  87                 m->errorOut(e, "GetMetaCommunityCommand", "GetMetaCommunityCommand");
  88                 exit(1);
  89         }
  90 }
  91 //**********************************************************************************************************************
  92 GetMetaCommunityCommand::GetMetaCommunityCommand(string option)  {
  93         try {
  94                 abort = false; calledHelp = false;
  95         allLines=true;
  96
  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;}
 100
 101                 else {
 102                         //valid paramters for this command
 103                         vector<string> myArray = setParameters();
 104
 105                         OptionParser parser(option);
 106                         map<string,string> parameters = parser.getParameters();
 107
 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;  }
 113                         }
 114
 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;
 122
 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 = "";          }
 126                         else {
 127                 string path;
 128                 it = parameters.find("shared");
 129                                 if(it != parameters.end()){
 130                                         path = m->hasPath(it->second);
 131                                         if (path == "") {       parameters["shared"] = inputDir + it->second;           }
 132                                 }
 133             }
 134
 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); }
 144
 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
 148                         }
 149
 150             string temp = validParameter.validFile(parameters, "minpartitions", false); if (temp == "not found"){       temp = "5";      }
 151                         m->mothurConvert(temp, minpartitions);
 152
 153             temp = validParameter.validFile(parameters, "maxpartitions", false);        if (temp == "not found"){       temp = "10";     }
 154                         m->mothurConvert(temp, maxpartitions);
 155
 156             temp = validParameter.validFile(parameters, "optimizegap", false);          if (temp == "not found"){       temp = "3";      }
 157                         m->mothurConvert(temp, optimizegap);
 158
 159             temp = validParameter.validFile(parameters, "processors", false);   if (temp == "not found"){       temp = m->getProcessors();      }
 160                         m->setProcessors(temp);
 161                         m->mothurConvert(temp, processors);
 162
 163             string groups = validParameter.validFile(parameters, "groups", false);
 164                         if (groups == "not found") { groups = ""; }
 165                         else { m->splitAtDash(groups, Groups); }
 166                         m->setGroups(Groups);
 167
 168             string label = validParameter.validFile(parameters, "label", false);
 169                         if (label == "not found") { label = ""; }
 170                         else {
 171                                 if(label != "all") {  m->splitAtDash(label, labels);  allLines = 0;  }
 172                                 else { allLines = 1;  }
 173                         }
 174                 }
 175
 176         }
 177         catch(exception& e) {
 178                 m->errorOut(e, "GetMetaCommunityCommand", "GetMetaCommunityCommand");
 179                 exit(1);
 180         }
 181 }
 182 //**********************************************************************************************************************
 183
 184 int GetMetaCommunityCommand::execute(){
 185         try {
 186
 187                 if (abort == true) { if (calledHelp) { return 0; }  return 2;   }
 188
 189         InputData input(sharedfile, "sharedfile");
 190         vector<SharedRAbundVector*> lookup = input.getSharedRAbundVectors();
 191         string lastLabel = lookup[0]->getLabel();
 192
 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;
 196
 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))) {
 199
 200             if (m->control_pressed) { for (int i = 0; i < lookup.size(); i++) {  delete lookup[i];  }  return 0; }
 201
 202             if(allLines == 1 || labels.count(lookup[0]->getLabel()) == 1){
 203
 204                 m->mothurOut(lookup[0]->getLabel()); m->mothurOutEndLine();
 205
 206                 createProcesses(lookup);
 207
 208                 processedLabels.insert(lookup[0]->getLabel());
 209                 userLabels.erase(lookup[0]->getLabel());
 210             }
 211
 212             if ((m->anyLabelsToProcess(lookup[0]->getLabel(), userLabels, "") == true) && (processedLabels.count(lastLabel) != 1)) {
 213                 string saveLabel = lookup[0]->getLabel();
 214
 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();
 218
 219                 createProcesses(lookup);
 220
 221                 processedLabels.insert(lookup[0]->getLabel());
 222                 userLabels.erase(lookup[0]->getLabel());
 223
 224                 //restore real lastlabel to save below
 225                 lookup[0]->setLabel(saveLabel);
 226             }
 227
 228             lastLabel = lookup[0]->getLabel();
 229             //prevent memory leak
 230             for (int i = 0; i < lookup.size(); i++) {  delete lookup[i]; lookup[i] = NULL; }
 231
 232             if (m->control_pressed) { return 0; }
 233
 234             //get next line to process
 235             lookup = input.getSharedRAbundVectors();
 236         }
 237
 238         if (m->control_pressed) {  return 0; }
 239
 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();
 247                 needToRun = true;
 248             }else {
 249                 m->mothurOut(". Please refer to " + lastLabel + "."); m->mothurOutEndLine();
 250             }
 251         }
 252
 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);
 257
 258             m->mothurOut(lookup[0]->getLabel()); m->mothurOutEndLine();
 259
 260             createProcesses(lookup);
 261
 262             for (int i = 0; i < lookup.size(); i++) {  delete lookup[i];  }
 263         }
 264
 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();
 270         return 0;
 271
 272     }
 273         catch(exception& e) {
 274                 m->errorOut(e, "GetMetaCommunityCommand", "execute");
 275                 exit(1);
 276         }
 277 }
 278 //**********************************************************************************************************************
 279 int GetMetaCommunityCommand::createProcesses(vector<SharedRAbundVector*>& thislookup){
 280         try {
 281
 282         #if defined (__APPLE__) || (__MACH__) || (linux) || (__linux) || (__linux__) || (__unix__) || (__unix)
 283         #else
 284         processors=1; //qFinderDMM not thread safe
 285         #endif
 286
 287         vector<int> processIDS;
 288                 int process = 1;
 289                 int num = 0;
 290         int minPartition = 0;
 291
 292                 //sanity check
 293                 if (maxpartitions < processors) { processors = maxpartitions; }
 294
 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);
 300
 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);
 308
 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++) {
 312             //cout << i << endl;
 313             int processToAssign = (i+1) % processors;
 314             if (processToAssign == 0) { processToAssign = processors; }
 315
 316             if (m->debug) { m->mothurOut("[DEBUG]: assigning " + toString(i) + " to process " + toString(processToAssign-1) + "\n"); }
 317             dividedPartitions[(processToAssign-1)].push_back(i);
 318
 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);
 324         }
 325
 326         for (int i = 0; i < processors; i++) { //read from everyone elses, just write to yours
 327             string tempDoneFile = toString(i) + ".done.temp";
 328             doneFlags.push_back(tempDoneFile);
 329             ofstream out;
 330             m->openOutputFile(tempDoneFile, out); //clear out
 331             out.close();
 332         }
 333
 334
 335 #if defined (__APPLE__) || (__MACH__) || (linux) || (__linux) || (__linux__) || (__unix__) || (__unix)
 336
 337                 //loop through and create all the processes you want
 338                 while (process != processors) {
 339                         int pid = fork();
 340
 341                         if (pid > 0) {
 342                                 processIDS.push_back(pid);  //create map from line number to pid so you can append files in correct order later
 343                                 process++;
 344                         }else if (pid == 0){
 345                 outputNames.clear();
 346                                 num = processDriver(thislookup, dividedPartitions[process], (outputFileName + toString(getpid())), rels[process], matrix[process], doneFlags, process);
 347
 348                 //pass numSeqs to parent
 349                                 ofstream out;
 350                                 string tempFile = toString(getpid()) + ".outputNames.temp";
 351                                 m->openOutputFile(tempFile, out);
 352                 out << num << endl;
 353                 out << outputNames.size() << endl;
 354                                 for (int i = 0; i < outputNames.size(); i++) { out << outputNames[i] << endl; }
 355                                 out.close();
 356
 357                                 exit(0);
 358                         }else {
 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); }
 361                                 exit(0);
 362                         }
 363                 }
 364
 365                 //do my part
 366                 minPartition = processDriver(thislookup, dividedPartitions[0], outputFileName, rels[0], matrix[0], doneFlags, 0);
 367
 368                 //force parent to wait until all the processes are done
 369                 for (int i=0;i<processIDS.size();i++) {
 370                         int temp = processIDS[i];
 371                         wait(&temp);
 372                 }
 373
 374         vector<string> tempOutputNames = outputNames;
 375         for (int i=0;i<processIDS.size();i++) {
 376             ifstream in;
 377                         string tempFile = toString(processIDS[i]) + ".outputNames.temp";
 378                         m->openInputFile(tempFile, in);
 379                         if (!in.eof()) {
 380                 int tempNum = 0;
 381                 in >> tempNum; m->gobble(in);
 382                 if (tempNum < minPartition) { minPartition = tempNum; }
 383                 in >> tempNum; m->gobble(in);
 384                 for (int i = 0; i < tempNum; i++) {
 385                     string tempName = "";
 386                     in >> tempName; m->gobble(in);
 387                     tempOutputNames.push_back(tempName);
 388                 }
 389             }
 390                         in.close(); m->mothurRemove(tempFile);
 391
 392             m->appendFilesWithoutHeaders(outputFileName + toString(processIDS[i]), outputFileName);
 393             m->mothurRemove(outputFileName + toString(processIDS[i]));
 394         }
 395
 396         if (processors > 1) {
 397             outputNames.clear();
 398             for (int i = 0; i < tempOutputNames.size(); i++) { //remove files if needed
 399                 string name = tempOutputNames[i];
 400                 vector<string> parts;
 401                 m->splitAtChar(name, parts, '.');
 402                 bool keep = true;
 403                 if (((parts[parts.size()-1] == "relabund") || (parts[parts.size()-1] == "posterior")) && (parts[parts.size()-2] == "mix")) {
 404                     string tempNum = parts[parts.size()-3];
 405                     int num;  m->mothurConvert(tempNum, num);
 406                     //if (num > minPartition) {
 407                      //   m->mothurRemove(tempOutputNames[i]);
 408                     //    keep = false; if (m->debug) { m->mothurOut("[DEBUG]: removing " + tempOutputNames[i] + ".\n"); }
 409                     //}
 410                 }
 411                 if (keep) { outputNames.push_back(tempOutputNames[i]); }
 412             }
 413
 414             //reorder fit file
 415             ifstream in;
 416             m->openInputFile(outputFileName, in);
 417             string headers = m->getline(in); m->gobble(in);
 418
 419             map<int, string> file;
 420             while (!in.eof()) {
 421                 string numString, line;
 422                 int num;
 423                 in >> numString; line = m->getline(in); m->gobble(in);
 424                 m->mothurConvert(numString, num);
 425                 file[num] = line;
 426             }
 427             in.close();
 428             ofstream out;
 429             m->openOutputFile(outputFileName, out);
 430             out << headers << endl;
 431             for (map<int, string>::iterator it = file.begin(); it != file.end(); it++) {
 432                 out << it->first << '\t' << it->second << endl;
 433                 if (m->debug) { m->mothurOut("[DEBUG]: printing: " + toString(it->first) + '\t' + it->second + ".\n"); }
 434             }
 435             out.close();
 436         }
 437
 438 #else
 439         /*
 440         vector<metaCommunityData*> pDataArray;
 441                 DWORD   dwThreadIdArray[processors-1];
 442                 HANDLE  hThreadArray[processors-1];
 443
 444                 //Create processor worker threads.
 445                 for( int i=1; i<processors; i++ ){
 446             //copy lookup
 447             //make copy of lookup so we don't get access violations
 448             vector<SharedRAbundVector*> newLookup;
 449             for (int k = 0; k < thislookup.size(); k++) {
 450                 SharedRAbundVector* temp = new SharedRAbundVector();
 451                 temp->setLabel(thislookup[k]->getLabel());
 452                 temp->setGroup(thislookup[k]->getGroup());
 453                 newLookup.push_back(temp);
 454             }
 455
 456             //for each bin
 457             for (int k = 0; k < thislookup[0]->getNumBins(); k++) {
 458                 if (m->control_pressed) { for (int j = 0; j < newLookup.size(); j++) {  delete newLookup[j];  } return 0; }
 459                 for (int j = 0; j < thislookup.size(); j++) { newLookup[j]->push_back(thislookup[j]->getAbundance(k), thislookup[j]->getGroup()); }
 460             }
 461
 462             processIDS.push_back(i);
 463
 464                         // Allocate memory for thread data.
 465                         metaCommunityData* tempMeta = new metaCommunityData(newLookup, m, dividedPartitions[i], outputFileName + toString(i), rels[i], matrix[i], minpartitions, maxpartitions, optimizegap);
 466                         pDataArray.push_back(tempMeta);
 467
 468                         hThreadArray[i] = CreateThread(NULL, 0, MyMetaCommunityThreadFunction, pDataArray[i], 0, &dwThreadIdArray[i]);
 469                 }
 470
 471         //do my part
 472                 minPartition = processDriver(thislookup, dividedPartitions[0], outputFileName, rels[0], matrix[0]);
 473
 474                 //Wait until all threads have terminated.
 475                 WaitForMultipleObjects(processors-1, hThreadArray, TRUE, INFINITE);
 476
 477                 //Close all thread handles and free memory allocations.
 478                 for(int i=0; i < pDataArray.size(); i++){
 479             if (pDataArray[i]->minPartition < minPartition) { minPartition = pDataArray[i]->minPartition; }
 480             for (int j = 0; j < pDataArray[i]->outputNames.size(); j++) {
 481                 outputNames.push_back(pDataArray[i]->outputNames[j]);
 482             }
 483             m->appendFilesWithoutHeaders(outputFileName + toString(processIDS[i]), outputFileName);
 484             m->mothurRemove(outputFileName + toString(processIDS[i]));
 485                         CloseHandle(hThreadArray[i]);
 486                         delete pDataArray[i];
 487                 } */
 488         //do my part
 489                 minPartition = processDriver(thislookup, dividedPartitions[0], outputFileName, rels[0], matrix[0], doneFlags, 0);
 490 #endif
 491         for (int i = 0; i < processors; i++) { //read from everyone elses, just write to yours
 492             string tempDoneFile = toString(i) + ".done.temp";
 493             m->mothurRemove(tempDoneFile);
 494         }
 495
 496         if (m->control_pressed) { return 0; }
 497
 498         if (m->debug) { m->mothurOut("[DEBUG]: minPartition = " + toString(minPartition) + "\n"); }
 499
 500         //run generate Summary function for smallest minPartition
 501         variables["[tag]"] = toString(minPartition);
 502         generateSummaryFile(minPartition, variables);
 503
 504         return 0;
 505
 506     }
 507         catch(exception& e) {
 508                 m->errorOut(e, "GetMetaCommunityCommand", "createProcesses");
 509                 exit(1);
 510         }
 511 }
 512 //**********************************************************************************************************************
 513 int GetMetaCommunityCommand::processDriver(vector<SharedRAbundVector*>& thislookup, vector<int>& parts, string outputFileName, vector<string> relabunds, vector<string> matrix, vector<string> doneFlags, int processID){
 514         try {
 515
 516         double minLaplace = 1e10;
 517         int minPartition = 0;
 518
 519                 ofstream fitData;
 520                 m->openOutputFile(outputFileName, fitData);
 521         fitData.setf(ios::fixed, ios::floatfield);
 522         fitData.setf(ios::showpoint);
 523         cout.setf(ios::fixed, ios::floatfield);
 524         cout.setf(ios::showpoint);
 525
 526         vector< vector<int> > sharedMatrix;
 527         for (int i = 0; i < thislookup.size(); i++) { sharedMatrix.push_back(thislookup[i]->getAbundances()); }
 528
 529         m->mothurOut("K\tNLE\t\tlogDet\tBIC\t\tAIC\t\tLaplace\n");
 530         fitData << "K\tNLE\tlogDet\tBIC\tAIC\tLaplace" << endl;
 531
 532         for(int i=0;i<parts.size();i++){
 533
 534             int numPartitions = parts[i];
 535
 536             if (m->debug) { m->mothurOut("[DEBUG]: running partition " + toString(numPartitions) + "\n"); }
 537
 538             if (m->control_pressed) { break; }
 539
 540             //check to see if anyone else is done
 541             for (int j = 0; j < doneFlags.size(); j++) {
 542                 if (!m->isBlank(doneFlags[j])) { //another process has finished
 543                     //are they done at a lower partition?
 544                     ifstream in;
 545                     m->openInputFile(doneFlags[j], in);
 546                     int tempNum;
 547                     in >> tempNum; in.close();
 548                     if (tempNum < numPartitions) { break; } //quit, because someone else has finished
 549                 }
 550             }
 551
 552             qFinderDMM findQ(sharedMatrix, numPartitions);
 553
 554             double laplace = findQ.getLaplace();
 555             m->mothurOut(toString(numPartitions) + '\t');
 556             cout << setprecision (2) << findQ.getNLL() << '\t' << findQ.getLogDet() << '\t';
 557             m->mothurOutJustToLog(toString(findQ.getNLL()) + '\t' + toString(findQ.getLogDet()) + '\t');
 558             cout << findQ.getBIC() << '\t' << findQ.getAIC() << '\t' << laplace;
 559             m->mothurOutJustToLog(toString(findQ.getBIC()) + '\t' + toString(findQ.getAIC()) + '\t' + toString(laplace));
 560
 561             fitData << numPartitions << '\t';
 562             fitData << setprecision (2) << findQ.getNLL() << '\t' << findQ.getLogDet() << '\t';
 563             fitData << findQ.getBIC() << '\t' << findQ.getAIC() << '\t' << laplace << endl;
 564
 565             if(laplace < minLaplace){
 566                 minPartition = numPartitions;
 567                 minLaplace = laplace;
 568                 m->mothurOut("***");
 569             }
 570             m->mothurOutEndLine();
 571
 572             string relabund = relabunds[i];
 573             outputNames.push_back(relabund); outputTypes["relabund"].push_back(relabund);
 574             string matrixName = matrix[i];
 575             outputNames.push_back(matrixName); outputTypes["matrix"].push_back(matrixName);
 576
 577             findQ.printZMatrix(matrixName, m->getGroups());
 578             findQ.printRelAbund(relabund, m->currentBinLabels);
 579
 580             if(optimizegap != -1 && (numPartitions - minPartition) >= optimizegap && numPartitions >= minpartitions){
 581                 string tempDoneFile = toString(processID) + ".done.temp";
 582                 ofstream outDone;
 583                 m->openOutputFile(tempDoneFile, outDone);
 584                 outDone << minPartition << endl;
 585                 outDone.close();
 586                 break;
 587             }
 588         }
 589         fitData.close();
 590
 591         //minPartition = 4;
 592
 593         if (m->control_pressed) { return 0; }
 594
 595         return minPartition;
 596     }
 597         catch(exception& e) {
 598                 m->errorOut(e, "GetMetaCommunityCommand", "processDriver");
 599                 exit(1);
 600         }
 601 }
 602 /**************************************************************************************************/
 603
 604 vector<double> GetMetaCommunityCommand::generateDesignFile(int numPartitions, map<string,string> variables){
 605     try {
 606         vector<double> piValues(numPartitions, 0);
 607
 608         ifstream postFile;
 609         variables["[tag]"] = toString(numPartitions);
 610         string input = getOutputFileName("matrix", variables);
 611         m->openInputFile(input, postFile);//((fileRoot + toString(numPartitions) + "mix.posterior").c_str()); //matrix file
 612         variables.erase("[tag]");
 613                 string outputFileName = getOutputFileName("design", variables);
 614         ofstream designFile;
 615         m->openOutputFile(outputFileName, designFile);
 616         outputNames.push_back(outputFileName); outputTypes["design"].push_back(outputFileName);
 617
 618
 619         vector<string> titles(numPartitions);
 620
 621         for(int i=0;i<numPartitions;i++){   postFile >> titles[i];  }
 622
 623         double posterior;
 624         string sampleName;
 625         int numSamples = 0;
 626
 627         while(postFile){
 628
 629             if (m->control_pressed) { break; }
 630
 631             double maxPosterior = 0.0000;
 632             int maxPartition = -1;
 633
 634             postFile >> sampleName;
 635
 636             for(int i=0;i<numPartitions;i++){
 637
 638                 postFile >> posterior;
 639                 if(posterior > maxPosterior){
 640                     maxPosterior = posterior;
 641                     maxPartition = i;
 642                 }
 643                 piValues[i] += posterior;
 644
 645             }
 646
 647             designFile << sampleName << '\t' << titles[maxPartition] << endl;
 648
 649             numSamples++;
 650             m->gobble(postFile);
 651         }
 652         for(int i=0;i<numPartitions;i++){
 653             piValues[i] /= (double)numSamples;
 654         }
 655
 656
 657         postFile.close();
 658         designFile.close();
 659
 660         return piValues;
 661     }
 662         catch(exception& e) {
 663                 m->errorOut(e, "GetMetaCommunityCommand", "generateDesignFile");
 664                 exit(1);
 665         }
 666 }
 667
 668 /**************************************************************************************************/
 669
 670 inline bool summaryFunction(summaryData i, summaryData j){ return i.difference > j.difference;   }
 671
 672 /**************************************************************************************************/
 673 int GetMetaCommunityCommand::generateSummaryFile(int numPartitions, map<string,string> v){
 674     try {
 675         vector<summaryData> summary;
 676
 677         vector<double> pMean(numPartitions, 0);
 678         vector<double> pLCI(numPartitions, 0);
 679         vector<double> pUCI(numPartitions, 0);
 680
 681         string name, header;
 682         double mean, lci, uci;
 683
 684
 685         vector<double> piValues = generateDesignFile(numPartitions, v);
 686
 687         ifstream referenceFile;
 688         map<string, string> variables;
 689         variables["[filename]"] = v["[filename]"];
 690         variables["[distance]"] = v["[distance]"];
 691         variables["[tag]"] = "1";
 692         string reference = getOutputFileName("relabund", variables);
 693         m->openInputFile(reference, referenceFile); //((fileRoot + label + ".1mix.relabund").c_str());
 694         variables["[tag]"] = toString(numPartitions);
 695         string partFile = getOutputFileName("relabund", variables);
 696         ifstream partitionFile;
 697         m->openInputFile(partFile, partitionFile); //((fileRoot + toString(numPartitions) + "mix.relabund").c_str());
 698
 699         header = m->getline(referenceFile);
 700         header = m->getline(partitionFile);
 701         stringstream head(header);
 702         string dummy, label;
 703         head >> dummy;
 704         vector<string> thetaValues(numPartitions, "");
 705         for(int i=0;i<numPartitions;i++){
 706             head >> label >> dummy >> dummy;
 707             thetaValues[i] = label.substr(label.find_last_of('_')+1);
 708         }
 709
 710
 711         vector<double> partitionDiff(numPartitions, 0.0000);
 712
 713         while(referenceFile){
 714
 715             if (m->control_pressed) { break; }
 716
 717             referenceFile >> name >> mean >> lci >> uci;
 718
 719             summaryData tempData;
 720             tempData.name = name;
 721             tempData.refMean = mean;
 722
 723             double difference = 0.0000;
 724
 725             partitionFile >> name;
 726             for(int j=0;j<numPartitions;j++){
 727                 partitionFile >> pMean[j] >> pLCI[j] >> pUCI[j];
 728                 difference += abs(mean - pMean[j]);
 729                 partitionDiff[j] += abs(mean - pMean[j]);;
 730             }
 731
 732             tempData.partMean = pMean;
 733             tempData.partLCI = pLCI;
 734             tempData.partUCI = pUCI;
 735             tempData.difference = difference;
 736             summary.push_back(tempData);
 737
 738             m->gobble(referenceFile);
 739             m->gobble(partitionFile);
 740         }
 741         referenceFile.close();
 742         partitionFile.close();
 743
 744         if (m->control_pressed) { return 0; }
 745
 746         int numOTUs = (int)summary.size();
 747
 748         sort(summary.begin(), summary.end(), summaryFunction);
 749
 750         variables.erase("[tag]");
 751                 string outputFileName = getOutputFileName("parameters", variables);
 752         outputNames.push_back(outputFileName); outputTypes["parameters"].push_back(outputFileName);
 753
 754         ofstream parameterFile;
 755         m->openOutputFile(outputFileName, parameterFile); //((fileRoot + "mix.parameters").c_str());
 756         parameterFile.setf(ios::fixed, ios::floatfield);
 757         parameterFile.setf(ios::showpoint);
 758
 759         double totalDifference =  0.0000;
 760         parameterFile << "Part\tDif2Ref_i\ttheta_i\tpi_i\n";
 761         for(int i=0;i<numPartitions;i++){
 762             if (m->control_pressed) { break; }
 763             parameterFile << i+1 << '\t' << setprecision(2) << partitionDiff[i] << '\t' << thetaValues[i] << '\t' << piValues[i] << endl;
 764             totalDifference += partitionDiff[i];
 765         }
 766         parameterFile.close();
 767
 768         if (m->control_pressed) { return 0; }
 769
 770         string summaryFileName = getOutputFileName("summary", variables);
 771         outputNames.push_back(summaryFileName); outputTypes["summary"].push_back(summaryFileName);
 772
 773         ofstream summaryFile;
 774         m->openOutputFile(summaryFileName, summaryFile); //((fileRoot + "mix.summary").c_str());
 775         summaryFile.setf(ios::fixed, ios::floatfield);
 776         summaryFile.setf(ios::showpoint);
 777
 778
 779         summaryFile << "OTU\tP0.mean";
 780         for(int i=0;i<numPartitions;i++){
 781             summaryFile << "\tP" << i+1 << ".mean\tP" << i+1 << ".lci\tP" << i+1 << ".uci";
 782         }
 783         summaryFile << "\tDifference\tCumFraction" << endl;
 784
 785         double cumDiff = 0.0000;
 786
 787         for(int i=0;i<numOTUs;i++){
 788             if (m->control_pressed) { break; }
 789             summaryFile << summary[i].name << setprecision(2) << '\t' << summary[i].refMean;
 790             for(int j=0;j<numPartitions;j++){
 791                 summaryFile  << '\t' << summary[i].partMean[j] << '\t' << summary[i].partLCI[j] << '\t' << summary[i].partUCI[j];
 792             }
 793
 794             cumDiff += summary[i].difference/totalDifference;
 795             summaryFile << '\t' << summary[i].difference << '\t' << cumDiff << endl;
 796         }
 797         summaryFile.close();
 798
 799         return 0;
 800
 801     }
 802         catch(exception& e) {
 803                 m->errorOut(e, "GetMetaCommunityCommand", "generateSummaryFile");
 804                 exit(1);
 805         }
 806
 807 }
 808 //**********************************************************************************************************************
 809