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 #include "communitytype.h"
  11 #include "kmeans.h"
  12 #include "validcalculator.h"
  13 #include "subsample.h"
  14
  15 //**********************************************************************************************************************
  16 vector<string> GetMetaCommunityCommand::setParameters(){
  17         try {
  18         CommandParameter pshared("shared", "InputTypes", "", "", "none", "none", "none","outputType",false,true); parameters.push_back(pshared);
  19         CommandParameter pgroups("groups", "String", "", "", "", "", "","",false,false); parameters.push_back(pgroups);
  20                 CommandParameter plabel("label", "String", "", "", "", "", "","",false,false); parameters.push_back(plabel);
  21         CommandParameter pcalc("calc", "Multiple", "sharedsobs-sharedchao-sharedace-jabund-sorabund-jclass-sorclass-jest-sorest-thetayc-thetan-kstest-sharednseqs-ochiai-anderberg-kulczynski-kulczynskicody-lennon-morisitahorn-braycurtis-whittaker-odum-canberra-structeuclidean-structchord-hellinger-manhattan-structpearson-soergel-spearman-structkulczynski-speciesprofile-hamming-structchi2-gower-memchi2-memchord-memeuclidean-mempearson-jsd-rjsd", "rjsd", "", "", "","",false,false,true); parameters.push_back(pcalc);
  22         CommandParameter psubsample("subsample", "String", "", "", "", "", "","",false,false); parameters.push_back(psubsample);
  23         CommandParameter piters("iters", "Number", "", "1000", "", "", "","",false,false); parameters.push_back(piters);
  24         CommandParameter pminpartitions("minpartitions", "Number", "", "5", "", "", "","",false,false,true); parameters.push_back(pminpartitions);
  25         CommandParameter pmaxpartitions("maxpartitions", "Number", "", "100", "", "", "","",false,false,true); parameters.push_back(pmaxpartitions);
  26         CommandParameter poptimizegap("optimizegap", "Number", "", "3", "", "", "","",false,false,true); parameters.push_back(poptimizegap);
  27         CommandParameter pprocessors("processors", "Number", "", "1", "", "", "","",false,false,true); parameters.push_back(pprocessors);
  28                 CommandParameter pinputdir("inputdir", "String", "", "", "", "", "","",false,false); parameters.push_back(pinputdir);
  29                 CommandParameter poutputdir("outputdir", "String", "", "", "", "", "","",false,false); parameters.push_back(poutputdir);
  30                 CommandParameter pmethod("method", "Multiple", "dmm-kmeans-pam", "dmm", "", "", "","",false,false,true); parameters.push_back(pmethod);
  31
  32                 vector<string> myArray;
  33                 for (int i = 0; i < parameters.size(); i++) {   myArray.push_back(parameters[i].name);          }
  34                 return myArray;
  35         }
  36         catch(exception& e) {
  37                 m->errorOut(e, "NewCommand", "setParameters");
  38                 exit(1);
  39         }
  40 }
  41 //**********************************************************************************************************************
  42 string GetMetaCommunityCommand::getHelpString(){
  43         try {
  44                 string helpString = "";
  45                 helpString += "The get.communitytype command parameters are shared, method, label, groups, minpartitions, maxpartitions, optimizegap and processors. The shared file is required. \n";
  46         helpString += "The label parameter is used to analyze specific labels in your input. labels are separated by dashes.\n";
  47                 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";
  48         helpString += "The method parameter allows to select the method you would like to use.  Options are dmm, kmeans and pam. Default=dmm.\n";
  49         helpString += "The calc parameter allows to select the calculator you would like to use to calculate the distance matrix used by the pam method. By default the jsd calculator is used.\n";
  50         helpString += "The iters parameter allows you to choose the number of times you would like to run the subsample while calculating the distance matirx for the pam method.\n";
  51         helpString += "The subsample parameter allows you to enter the size pergroup of the sample or you can set subsample=T and mothur will use the size of your smallest group while calculating the distance matrix for the pam method.\n";
  52                 helpString += "The minpartitions parameter is used to .... Default=5.\n";
  53         helpString += "The maxpartitions parameter is used to .... Default=10.\n";
  54         helpString += "The optimizegap parameter is used to .... Default=3.\n";
  55         helpString += "The processors parameter allows you to specify number of processors to use.  The default is 1.\n";
  56                 helpString += "The get.communitytype command should be in the following format: get.communitytype(shared=yourSharedFile).\n";
  57                 return helpString;
  58         }
  59         catch(exception& e) {
  60                 m->errorOut(e, "GetMetaCommunityCommand", "getHelpString");
  61                 exit(1);
  62         }
  63 }
  64 //**********************************************************************************************************************
  65 string GetMetaCommunityCommand::getOutputPattern(string type) {
  66     try {
  67         string pattern = "";
  68
  69         if (type == "fit")              {  pattern = "[filename],[distance],mix.fit"; }
  70         else if (type == "relabund")    {  pattern = "[filename],[distance],[tag],mix.relabund"; }
  71         else if (type == "design")      {  pattern = "[filename],[distance],mix.design"; }
  72         else if (type == "matrix")      {  pattern = "[filename],[distance],[tag],mix.posterior"; }
  73         else if (type == "parameters")  {  pattern = "[filename],[distance],mix.parameters"; }
  74         else if (type == "summary")  {  pattern = "[filename],[distance],mix.summary"; }
  75         else { m->mothurOut("[ERROR]: No definition for type " + type + " output pattern.\n"); m->control_pressed = true;  }
  76
  77         return pattern;
  78     }
  79     catch(exception& e) {
  80         m->errorOut(e, "GetMetaCommunityCommand", "getOutputPattern");
  81         exit(1);
  82     }
  83 }
  84 //**********************************************************************************************************************
  85 GetMetaCommunityCommand::GetMetaCommunityCommand(){
  86         try {
  87                 abort = true; calledHelp = true;
  88                 setParameters();
  89         vector<string> tempOutNames;
  90                 outputTypes["fit"] = tempOutNames;
  91         outputTypes["relabund"] = tempOutNames;
  92         outputTypes["matrix"] = tempOutNames;
  93         outputTypes["design"] = tempOutNames;
  94         outputTypes["parameters"] = tempOutNames;
  95         outputTypes["summary"] = tempOutNames;
  96         }
  97         catch(exception& e) {
  98                 m->errorOut(e, "GetMetaCommunityCommand", "GetMetaCommunityCommand");
  99                 exit(1);
 100         }
 101 }
 102 //**********************************************************************************************************************
 103 GetMetaCommunityCommand::GetMetaCommunityCommand(string option)  {
 104         try {
 105                 abort = false; calledHelp = false;
 106         allLines=true;
 107
 108                 //allow user to run help
 109                 if(option == "help") { help(); abort = true; calledHelp = true; }
 110                 else if(option == "citation") { citation(); abort = true; calledHelp = true;}
 111
 112                 else {
 113                         //valid paramters for this command
 114                         vector<string> myArray = setParameters();
 115
 116                         OptionParser parser(option);
 117                         map<string,string> parameters = parser.getParameters();
 118
 119                         ValidParameters validParameter;
 120                         map<string,string>::iterator it;
 121                         //check to make sure all parameters are valid for command
 122                         for (it = parameters.begin(); it != parameters.end(); it++) {
 123                                 if (validParameter.isValidParameter(it->first, myArray, it->second) != true) {  abort = true;  }
 124                         }
 125
 126             vector<string> tempOutNames;
 127             outputTypes["fit"] = tempOutNames;
 128             outputTypes["relabund"] = tempOutNames;
 129             outputTypes["matrix"] = tempOutNames;
 130             outputTypes["design"] = tempOutNames;
 131             outputTypes["parameters"] = tempOutNames;
 132                         outputTypes["summary"] = tempOutNames;
 133
 134                         //if the user changes the input directory command factory will send this info to us in the output parameter
 135                         string inputDir = validParameter.validFile(parameters, "inputdir", false);
 136                         if (inputDir == "not found"){   inputDir = "";          }
 137                         else {
 138                 string path;
 139                 it = parameters.find("shared");
 140                                 if(it != parameters.end()){
 141                                         path = m->hasPath(it->second);
 142                                         if (path == "") {       parameters["shared"] = inputDir + it->second;           }
 143                                 }
 144             }
 145
 146             //get shared file, it is required
 147                         sharedfile = validParameter.validFile(parameters, "shared", true);
 148                         if (sharedfile == "not open") { sharedfile = ""; abort = true; }
 149                         else if (sharedfile == "not found") {
 150                                 //if there is a current shared file, use it
 151                                 sharedfile = m->getSharedFile();
 152                                 if (sharedfile != "") { m->mothurOut("Using " + sharedfile + " as input file for the shared parameter."); m->mothurOutEndLine(); }
 153                                 else {  m->mothurOut("You have no current sharedfile and the shared parameter is required."); m->mothurOutEndLine(); abort = true; }
 154                         }else { m->setSharedFile(sharedfile); }
 155
 156             //if the user changes the output directory command factory will send this info to us in the output parameter
 157                         outputDir = validParameter.validFile(parameters, "outputdir", false);           if (outputDir == "not found"){
 158                                 outputDir = m->hasPath(sharedfile); //if user entered a file with a path then preserve it
 159                         }
 160
 161             string temp = validParameter.validFile(parameters, "minpartitions", false); if (temp == "not found"){       temp = "5";      }
 162                         m->mothurConvert(temp, minpartitions);
 163
 164             temp = validParameter.validFile(parameters, "maxpartitions", false);        if (temp == "not found"){       temp = "10";     }
 165                         m->mothurConvert(temp, maxpartitions);
 166
 167             temp = validParameter.validFile(parameters, "optimizegap", false);          if (temp == "not found"){       temp = "3";      }
 168                         m->mothurConvert(temp, optimizegap);
 169
 170             temp = validParameter.validFile(parameters, "processors", false);   if (temp == "not found"){       temp = m->getProcessors();      }
 171                         m->setProcessors(temp);
 172                         m->mothurConvert(temp, processors);
 173
 174             string groups = validParameter.validFile(parameters, "groups", false);
 175                         if (groups == "not found") { groups = ""; }
 176                         else { m->splitAtDash(groups, Groups); }
 177                         m->setGroups(Groups);
 178
 179             string label = validParameter.validFile(parameters, "label", false);
 180                         if (label == "not found") { label = ""; }
 181                         else {
 182                                 if(label != "all") {  m->splitAtDash(label, labels);  allLines = 0;  }
 183                                 else { allLines = 1;  }
 184                         }
 185
 186             method = validParameter.validFile(parameters, "method", false);
 187                         if (method == "not found") { method = "dmm"; }
 188
 189                         if ((method == "dmm") || (method == "kmeans") || (method == "pam")) { }
 190                         else { m->mothurOut("[ERROR]: " + method + " is not a valid method.  Valid algorithms are dmm, kmeans and pam."); m->mothurOutEndLine(); abort = true; }
 191
 192             calc = validParameter.validFile(parameters, "calc", false);
 193                         if (calc == "not found") { calc = "rjsd";  }
 194                         else {
 195                 if (calc == "default")  {  calc = "rjsd";  }
 196                         }
 197                         m->splitAtDash(calc, Estimators);
 198                         if (m->inUsersGroups("citation", Estimators)) {
 199                                 ValidCalculators validCalc; validCalc.printCitations(Estimators);
 200                                 //remove citation from list of calcs
 201                                 for (int i = 0; i < Estimators.size(); i++) { if (Estimators[i] == "citation") {  Estimators.erase(Estimators.begin()+i); break; } }
 202                         }
 203             if (Estimators.size() != 1) { abort = true; m->mothurOut("[ERROR]: only one calculator is allowed.\n"); }
 204
 205             temp = validParameter.validFile(parameters, "iters", false);                        if (temp == "not found") { temp = "1000"; }
 206                         m->mothurConvert(temp, iters);
 207
 208             temp = validParameter.validFile(parameters, "subsample", false);            if (temp == "not found") { temp = "F"; }
 209                         if (m->isNumeric1(temp)) { m->mothurConvert(temp, subsampleSize); subsample = true; }
 210             else {
 211                 if (m->isTrue(temp)) { subsample = true; subsampleSize = -1; }  //we will set it to smallest group later
 212                 else { subsample = false; }
 213             }
 214
 215             if (subsample == false) { iters = 0; }
 216                 }
 217
 218         }
 219         catch(exception& e) {
 220                 m->errorOut(e, "GetMetaCommunityCommand", "GetMetaCommunityCommand");
 221                 exit(1);
 222         }
 223 }
 224 //**********************************************************************************************************************
 225
 226 int GetMetaCommunityCommand::execute(){
 227         try {
 228
 229                 if (abort == true) { if (calledHelp) { return 0; }  return 2;   }
 230
 231         InputData input(sharedfile, "sharedfile");
 232         vector<SharedRAbundVector*> lookup = input.getSharedRAbundVectors();
 233         string lastLabel = lookup[0]->getLabel();
 234
 235         //if the users enters label "0.06" and there is no "0.06" in their file use the next lowest label.
 236         set<string> processedLabels;
 237         set<string> userLabels = labels;
 238
 239         if (subsample) {
 240             if (subsampleSize == -1) { //user has not set size, set size = smallest samples size
 241                 subsampleSize = lookup[0]->getNumSeqs();
 242                 for (int i = 1; i < lookup.size(); i++) {
 243                     int thisSize = lookup[i]->getNumSeqs();
 244
 245                     if (thisSize < subsampleSize) {     subsampleSize = thisSize;       }
 246                 }
 247             }else {
 248                 m->clearGroups();
 249                 Groups.clear();
 250                 vector<SharedRAbundVector*> temp;
 251                 for (int i = 0; i < lookup.size(); i++) {
 252                     if (lookup[i]->getNumSeqs() < subsampleSize) {
 253                         m->mothurOut(lookup[i]->getGroup() + " contains " + toString(lookup[i]->getNumSeqs()) + ". Eliminating."); m->mothurOutEndLine();
 254                         delete lookup[i];
 255                     }else {
 256                         Groups.push_back(lookup[i]->getGroup());
 257                         temp.push_back(lookup[i]);
 258                     }
 259                 }
 260                 lookup = temp;
 261                 m->setGroups(Groups);
 262             }
 263
 264             if (lookup.size() < 2) { m->mothurOut("You have not provided enough valid groups.  I cannot run the command."); m->mothurOutEndLine(); m->control_pressed = true;  return 0; }
 265         }
 266
 267
 268         //as long as you are not at the end of the file or done wih the lines you want
 269         while((lookup[0] != NULL) && ((allLines == 1) || (userLabels.size() != 0))) {
 270
 271             if (m->control_pressed) { for (int i = 0; i < lookup.size(); i++) {  delete lookup[i];  }  return 0; }
 272
 273             if(allLines == 1 || labels.count(lookup[0]->getLabel()) == 1){
 274
 275                 m->mothurOut(lookup[0]->getLabel()); m->mothurOutEndLine();
 276
 277                 createProcesses(lookup);
 278
 279                 processedLabels.insert(lookup[0]->getLabel());
 280                 userLabels.erase(lookup[0]->getLabel());
 281             }
 282
 283             if ((m->anyLabelsToProcess(lookup[0]->getLabel(), userLabels, "") == true) && (processedLabels.count(lastLabel) != 1)) {
 284                 string saveLabel = lookup[0]->getLabel();
 285
 286                 for (int i = 0; i < lookup.size(); i++) {  delete lookup[i];  }
 287                 lookup = input.getSharedRAbundVectors(lastLabel);
 288                 m->mothurOut(lookup[0]->getLabel()); m->mothurOutEndLine();
 289
 290                 createProcesses(lookup);
 291
 292                 processedLabels.insert(lookup[0]->getLabel());
 293                 userLabels.erase(lookup[0]->getLabel());
 294
 295                 //restore real lastlabel to save below
 296                 lookup[0]->setLabel(saveLabel);
 297             }
 298
 299             lastLabel = lookup[0]->getLabel();
 300             //prevent memory leak
 301             for (int i = 0; i < lookup.size(); i++) {  delete lookup[i]; lookup[i] = NULL; }
 302
 303             if (m->control_pressed) { return 0; }
 304
 305             //get next line to process
 306             lookup = input.getSharedRAbundVectors();
 307         }
 308
 309         if (m->control_pressed) {  return 0; }
 310
 311         //output error messages about any remaining user labels
 312         set<string>::iterator it;
 313         bool needToRun = false;
 314         for (it = userLabels.begin(); it != userLabels.end(); it++) {
 315             m->mothurOut("Your file does not include the label " + *it);
 316             if (processedLabels.count(lastLabel) != 1) {
 317                 m->mothurOut(". I will use " + lastLabel + "."); m->mothurOutEndLine();
 318                 needToRun = true;
 319             }else {
 320                 m->mothurOut(". Please refer to " + lastLabel + "."); m->mothurOutEndLine();
 321             }
 322         }
 323
 324         //run last label if you need to
 325         if (needToRun == true)  {
 326             for (int i = 0; i < lookup.size(); i++) { if (lookup[i] != NULL) { delete lookup[i]; } }
 327             lookup = input.getSharedRAbundVectors(lastLabel);
 328
 329             m->mothurOut(lookup[0]->getLabel()); m->mothurOutEndLine();
 330
 331             createProcesses(lookup);
 332
 333             for (int i = 0; i < lookup.size(); i++) {  delete lookup[i];  }
 334         }
 335
 336         //output files created by command
 337                 m->mothurOutEndLine();
 338                 m->mothurOut("Output File Names: "); m->mothurOutEndLine();
 339                 for (int i = 0; i < outputNames.size(); i++) {  m->mothurOut(outputNames[i]); m->mothurOutEndLine();    }
 340                 m->mothurOutEndLine();
 341         return 0;
 342
 343     }
 344         catch(exception& e) {
 345                 m->errorOut(e, "GetMetaCommunityCommand", "execute");
 346                 exit(1);
 347         }
 348 }
 349 //**********************************************************************************************************************
 350 int GetMetaCommunityCommand::createProcesses(vector<SharedRAbundVector*>& thislookup){
 351         try {
 352
 353         #if defined (__APPLE__) || (__MACH__) || (linux) || (__linux) || (__linux__) || (__unix__) || (__unix)
 354         #else
 355         processors=1; //qFinderDMM not thread safe
 356         #endif
 357
 358         vector<int> processIDS;
 359                 int process = 1;
 360                 int num = 0;
 361         int minPartition = 0;
 362
 363                 //sanity check
 364                 if (maxpartitions < processors) { processors = maxpartitions; }
 365
 366         map<string, string> variables;
 367         variables["[filename]"] = outputDir + m->getRootName(m->getSimpleName(sharedfile));
 368         variables["[distance]"] = thislookup[0]->getLabel();
 369                 string outputFileName = getOutputFileName("fit", variables);
 370         outputNames.push_back(outputFileName); outputTypes["fit"].push_back(outputFileName);
 371
 372                 //divide the partitions between the processors
 373                 vector< vector<int> > dividedPartitions;
 374         vector< vector<string> > rels, matrix;
 375         vector<string> doneFlags;
 376         dividedPartitions.resize(processors);
 377         rels.resize(processors);
 378         matrix.resize(processors);
 379
 380         //for each file group figure out which process will complete it
 381         //want to divide the load intelligently so the big files are spread between processes
 382         for (int i=1; i<=maxpartitions; i++) {
 383             //cout << i << endl;
 384             int processToAssign = (i+1) % processors;
 385             if (processToAssign == 0) { processToAssign = processors; }
 386
 387             if (m->debug) { m->mothurOut("[DEBUG]: assigning " + toString(i) + " to process " + toString(processToAssign-1) + "\n"); }
 388             dividedPartitions[(processToAssign-1)].push_back(i);
 389
 390             variables["[tag]"] = toString(i);
 391             string relName = getOutputFileName("relabund", variables);
 392             string mName = getOutputFileName("matrix", variables);
 393             rels[(processToAssign-1)].push_back(relName);
 394             matrix[(processToAssign-1)].push_back(mName);
 395         }
 396
 397         for (int i = 0; i < processors; i++) { //read from everyone elses, just write to yours
 398             string tempDoneFile = m->getRootName(m->getSimpleName(sharedfile)) + toString(i) + ".done.temp";
 399             doneFlags.push_back(tempDoneFile);
 400             ofstream out;
 401             m->openOutputFile(tempDoneFile, out); //clear out
 402             out.close();
 403         }
 404
 405
 406 #if defined (__APPLE__) || (__MACH__) || (linux) || (__linux) || (__linux__) || (__unix__) || (__unix)
 407
 408                 //loop through and create all the processes you want
 409                 while (process != processors) {
 410                         int pid = fork();
 411
 412                         if (pid > 0) {
 413                                 processIDS.push_back(pid);  //create map from line number to pid so you can append files in correct order later
 414                                 process++;
 415                         }else if (pid == 0){
 416                 outputNames.clear();
 417                                 num = processDriver(thislookup, dividedPartitions[process], (outputFileName + toString(getpid())), rels[process], matrix[process], doneFlags, process);
 418
 419                 //pass numSeqs to parent
 420                                 ofstream out;
 421                                 string tempFile = toString(getpid()) + ".outputNames.temp";
 422                                 m->openOutputFile(tempFile, out);
 423                 out << num << endl;
 424                 out << outputNames.size() << endl;
 425                                 for (int i = 0; i < outputNames.size(); i++) { out << outputNames[i] << endl; }
 426                                 out.close();
 427
 428                                 exit(0);
 429                         }else {
 430                                 m->mothurOut("[ERROR]: unable to spawn the necessary processes."); m->mothurOutEndLine();
 431                                 for (int i = 0; i < processIDS.size(); i++) { kill (processIDS[i], SIGINT); }
 432                                 exit(0);
 433                         }
 434                 }
 435
 436                 //do my part
 437         if (method == "dmm") {  m->mothurOut("K\tNLE\t\tlogDet\tBIC\t\tAIC\t\tLaplace\n");  }
 438         else {
 439             m->mothurOut("K\tCH\t");
 440             for (int i = 0; i < thislookup.size(); i++) {  m->mothurOut(thislookup[i]->getGroup() + '\t'); }
 441             m->mothurOut("\n");
 442         }
 443                 minPartition = processDriver(thislookup, dividedPartitions[0], outputFileName, rels[0], matrix[0], doneFlags, 0);
 444
 445                 //force parent to wait until all the processes are done
 446                 for (int i=0;i<processIDS.size();i++) {
 447                         int temp = processIDS[i];
 448                         wait(&temp);
 449                 }
 450
 451         vector<string> tempOutputNames = outputNames;
 452         for (int i=0;i<processIDS.size();i++) {
 453             ifstream in;
 454                         string tempFile = toString(processIDS[i]) + ".outputNames.temp";
 455                         m->openInputFile(tempFile, in);
 456                         if (!in.eof()) {
 457                 int tempNum = 0;
 458                 in >> tempNum; m->gobble(in);
 459                 if (tempNum < minPartition) { minPartition = tempNum; }
 460                 in >> tempNum; m->gobble(in);
 461                 for (int i = 0; i < tempNum; i++) {
 462                     string tempName = "";
 463                     in >> tempName; m->gobble(in);
 464                     tempOutputNames.push_back(tempName);
 465                 }
 466             }
 467                         in.close(); m->mothurRemove(tempFile);
 468
 469             m->appendFilesWithoutHeaders(outputFileName + toString(processIDS[i]), outputFileName);
 470             m->mothurRemove(outputFileName + toString(processIDS[i]));
 471         }
 472
 473         if (processors > 1) {
 474             outputNames.clear();
 475             for (int i = 0; i < tempOutputNames.size(); i++) { //remove files if needed
 476                 string name = tempOutputNames[i];
 477                 vector<string> parts;
 478                 m->splitAtChar(name, parts, '.');
 479                 bool keep = true;
 480                 if (((parts[parts.size()-1] == "relabund") || (parts[parts.size()-1] == "posterior")) && (parts[parts.size()-2] == "mix")) {
 481                     string tempNum = parts[parts.size()-3];
 482                     int num;  m->mothurConvert(tempNum, num);
 483                     //if (num > minPartition) {
 484                      //   m->mothurRemove(tempOutputNames[i]);
 485                     //    keep = false; if (m->debug) { m->mothurOut("[DEBUG]: removing " + tempOutputNames[i] + ".\n"); }
 486                     //}
 487                 }
 488                 if (keep) { outputNames.push_back(tempOutputNames[i]); }
 489             }
 490
 491             //reorder fit file
 492             ifstream in;
 493             m->openInputFile(outputFileName, in);
 494             string headers = m->getline(in); m->gobble(in);
 495
 496             map<int, string> file;
 497             while (!in.eof()) {
 498                 string numString, line;
 499                 int num;
 500                 in >> numString; line = m->getline(in); m->gobble(in);
 501                 m->mothurConvert(numString, num);
 502                 file[num] = line;
 503             }
 504             in.close();
 505             ofstream out;
 506             m->openOutputFile(outputFileName, out);
 507             out << headers << endl;
 508             for (map<int, string>::iterator it = file.begin(); it != file.end(); it++) {
 509                 out << it->first << '\t' << it->second << endl;
 510                 if (m->debug) { m->mothurOut("[DEBUG]: printing: " + toString(it->first) + '\t' + it->second + ".\n"); }
 511             }
 512             out.close();
 513         }
 514
 515 #else
 516         m->mothurOut("K\tNLE\t\tlogDet\tBIC\t\tAIC\t\tLaplace\n");
 517                 minPartition = processDriver(thislookup, dividedPartitions[0], outputFileName, rels[0], matrix[0], doneFlags, 0);
 518 #endif
 519         for (int i = 0; i < processors; i++) { //read from everyone elses, just write to yours
 520             string tempDoneFile = m->getRootName(m->getSimpleName(sharedfile)) + toString(i) + ".done.temp";
 521             m->mothurRemove(tempDoneFile);
 522         }
 523
 524         if (m->control_pressed) { return 0; }
 525
 526         if (m->debug) { m->mothurOut("[DEBUG]: minPartition = " + toString(minPartition) + "\n"); }
 527
 528         //run generate Summary function for smallest minPartition
 529         variables["[tag]"] = toString(minPartition);
 530         vector<double> piValues = generateDesignFile(minPartition, variables);
 531         if (method == "dmm") {  generateSummaryFile(minPartition, variables, piValues); } //pam doesn't make a relabund file
 532
 533         return 0;
 534
 535     }
 536         catch(exception& e) {
 537                 m->errorOut(e, "GetMetaCommunityCommand", "createProcesses");
 538                 exit(1);
 539         }
 540 }
 541 //**********************************************************************************************************************
 542 int GetMetaCommunityCommand::processDriver(vector<SharedRAbundVector*>& thislookup, vector<int>& parts, string outputFileName, vector<string> relabunds, vector<string> matrix, vector<string> doneFlags, int processID){
 543         try {
 544
 545         double minLaplace = 1e10;
 546         int minPartition = 1;
 547         vector<double> minSilhouettes; minSilhouettes.resize(thislookup.size(), 0);
 548
 549                 ofstream fitData, silData;
 550         if (method == "dmm") {
 551             m->openOutputFile(outputFileName, fitData);
 552             fitData.setf(ios::fixed, ios::floatfield);
 553             fitData.setf(ios::showpoint);
 554             fitData << "K\tNLE\tlogDet\tBIC\tAIC\tLaplace" << endl;
 555         }else if((method == "pam") || (method == "kmeans")) { //because ch is looking of maximal value
 556             minLaplace = 0;
 557             m->openOutputFile(outputFileName, silData);
 558             silData.setf(ios::fixed, ios::floatfield);
 559             silData.setf(ios::showpoint);
 560             silData << "K\tCH\t";
 561             for (int i = 0; i < thislookup.size(); i++) { silData << thislookup[i]->getGroup() << '\t';  }
 562             silData << endl;
 563         }
 564
 565         cout.setf(ios::fixed, ios::floatfield);
 566         cout.setf(ios::showpoint);
 567
 568         vector< vector<int> > sharedMatrix;
 569         vector<string> thisGroups;
 570         for (int i = 0; i < thislookup.size(); i++) { sharedMatrix.push_back(thislookup[i]->getAbundances()); thisGroups.push_back(thislookup[i]->getGroup()); }
 571
 572         vector< vector<double> > dists; //do we want to output this matrix??
 573         if ((method == "pam") || (method == "kmeans")) {  dists = generateDistanceMatrix(thislookup);  }
 574
 575         if (m->debug) {
 576             m->mothurOut("[DEBUG]: dists = \n");
 577             for (int i = 0; i < dists.size(); i++) {
 578                 if (m->control_pressed) { break; }
 579                 m->mothurOut("[DEBUG]: i = " + toString(i) + '\t');
 580                 for (int j = 0; j < i; j++) { m->mothurOut(toString(dists[i][j]) +"\t"); }
 581                 m->mothurOut("\n");
 582             }
 583         }
 584
 585         for(int i=0;i<parts.size();i++){
 586
 587             int numPartitions = parts[i];
 588
 589             if (m->debug) { m->mothurOut("[DEBUG]: running partition " + toString(numPartitions) + "\n"); }
 590
 591             if (m->control_pressed) { break; }
 592
 593             //check to see if anyone else is done
 594             for (int j = 0; j < doneFlags.size(); j++) {
 595                 if (!m->isBlank(doneFlags[j])) { //another process has finished
 596                     //are they done at a lower partition?
 597                     ifstream in;
 598                     m->openInputFile(doneFlags[j], in);
 599                     int tempNum;
 600                     in >> tempNum; in.close();
 601                     if (tempNum < numPartitions) { break; } //quit, because someone else has finished
 602                 }
 603             }
 604
 605             CommunityTypeFinder* finder = NULL;
 606             if (method == "dmm")            {   finder = new qFinderDMM(sharedMatrix, numPartitions);   }
 607             else if (method == "kmeans")    {   finder = new KMeans(sharedMatrix, numPartitions);       }
 608             else if (method == "pam")       {   finder = new Pam(sharedMatrix, dists, numPartitions);                 }
 609             else {
 610                 if (i == 0) {  m->mothurOut(method + " is not a valid method option. I will run the command using dmm.\n"); }
 611                 finder = new qFinderDMM(sharedMatrix, numPartitions);
 612             }
 613
 614             double chi; vector<double> silhouettes;
 615             if (method == "dmm") {
 616                 double laplace = finder->getLaplace();
 617                 if(laplace < minLaplace){
 618                     minPartition = numPartitions;
 619                     minLaplace = laplace;
 620                 }
 621             }else {
 622                 chi = finder->calcCHIndex(dists);
 623                 silhouettes = finder->calcSilhouettes(dists);
 624                 if (chi > minLaplace) { //save partition with maximum ch index score
 625                     minPartition = numPartitions;
 626                     minLaplace = chi;
 627                     minSilhouettes = silhouettes;
 628                 }
 629             }
 630             string relabund = relabunds[i];
 631             string matrixName = matrix[i];
 632             outputNames.push_back(matrixName); outputTypes["matrix"].push_back(matrixName);
 633
 634             finder->printZMatrix(matrixName, thisGroups);
 635
 636             if (method == "dmm") {
 637                 finder->printFitData(cout, minLaplace);
 638                 finder->printFitData(fitData);
 639                 finder->printRelAbund(relabund, m->currentSharedBinLabels);
 640                 outputNames.push_back(relabund); outputTypes["relabund"].push_back(relabund);
 641             }else if ((method == "pam") || (method == "kmeans")) { //print silouettes and ch values
 642                 finder->printSilData(cout, chi, silhouettes);
 643                 finder->printSilData(silData, chi, silhouettes);
 644                 if (method == "kmeans") {
 645                     finder->printRelAbund(relabund, m->currentSharedBinLabels);
 646                     outputNames.push_back(relabund); outputTypes["relabund"].push_back(relabund);
 647                 }
 648             }
 649             delete finder;
 650
 651             if(optimizegap != -1 && (numPartitions - minPartition) >= optimizegap && numPartitions >= minpartitions){
 652                 string tempDoneFile = m->getRootName(m->getSimpleName(sharedfile)) + toString(processID) + ".done.temp";
 653                 ofstream outDone;
 654                 m->openOutputFile(tempDoneFile, outDone);
 655                 outDone << minPartition << endl;
 656                 outDone.close();
 657                 break;
 658             }
 659         }
 660         if (method == "dmm") { fitData.close(); }
 661
 662         if (m->control_pressed) { return 0; }
 663
 664         return minPartition;
 665     }
 666         catch(exception& e) {
 667                 m->errorOut(e, "GetMetaCommunityCommand", "processDriver");
 668                 exit(1);
 669         }
 670 }
 671 /**************************************************************************************************/
 672
 673 vector<double> GetMetaCommunityCommand::generateDesignFile(int numPartitions, map<string,string> variables){
 674     try {
 675         vector<double> piValues(numPartitions, 0);
 676
 677         ifstream postFile;
 678         variables["[tag]"] = toString(numPartitions);
 679         string input = getOutputFileName("matrix", variables);
 680         m->openInputFile(input, postFile);//((fileRoot + toString(numPartitions) + "mix.posterior").c_str()); //matrix file
 681         variables.erase("[tag]");
 682                 string outputFileName = getOutputFileName("design", variables);
 683         ofstream designFile;
 684         m->openOutputFile(outputFileName, designFile);
 685         outputNames.push_back(outputFileName); outputTypes["design"].push_back(outputFileName);
 686
 687
 688         vector<string> titles(numPartitions);
 689
 690         for(int i=0;i<numPartitions;i++){   postFile >> titles[i];  }
 691
 692         double posterior;
 693         string sampleName;
 694         int numSamples = 0;
 695
 696         while(postFile){
 697
 698             if (m->control_pressed) { break; }
 699
 700             double maxPosterior = 0.0000;
 701             int maxPartition = -1;
 702
 703             postFile >> sampleName;
 704
 705             for(int i=0;i<numPartitions;i++){
 706
 707                 postFile >> posterior;
 708                 if(posterior > maxPosterior){
 709                     maxPosterior = posterior;
 710                     maxPartition = i;
 711                 }
 712                 piValues[i] += posterior;
 713
 714             }
 715
 716             designFile << sampleName << '\t' << titles[maxPartition] << endl;
 717
 718             numSamples++;
 719             m->gobble(postFile);
 720         }
 721         for(int i=0;i<numPartitions;i++){
 722             piValues[i] /= (double)numSamples;
 723         }
 724
 725
 726         postFile.close();
 727         designFile.close();
 728
 729         return piValues;
 730     }
 731         catch(exception& e) {
 732                 m->errorOut(e, "GetMetaCommunityCommand", "generateDesignFile");
 733                 exit(1);
 734         }
 735 }
 736
 737 /**************************************************************************************************/
 738
 739 inline bool summaryFunction(summaryData i, summaryData j){ return i.difference > j.difference;   }
 740
 741 /**************************************************************************************************/
 742 int GetMetaCommunityCommand::generateSummaryFile(int numPartitions, map<string,string> v, vector<double> piValues){
 743     try {
 744         vector<summaryData> summary;
 745
 746         vector<double> pMean(numPartitions, 0);
 747         vector<double> pLCI(numPartitions, 0);
 748         vector<double> pUCI(numPartitions, 0);
 749
 750         string name, header;
 751         double mean, lci, uci;
 752
 753         ifstream referenceFile;
 754         map<string, string> variables;
 755         variables["[filename]"] = v["[filename]"];
 756         variables["[distance]"] = v["[distance]"];
 757         variables["[tag]"] = "1";
 758         string reference = getOutputFileName("relabund", variables);
 759         m->openInputFile(reference, referenceFile); //((fileRoot + label + ".1mix.relabund").c_str());
 760         variables["[tag]"] = toString(numPartitions);
 761         string partFile = getOutputFileName("relabund", variables);
 762         ifstream partitionFile;
 763         m->openInputFile(partFile, partitionFile); //((fileRoot + toString(numPartitions) + "mix.relabund").c_str());
 764
 765         header = m->getline(referenceFile);
 766         header = m->getline(partitionFile);
 767         stringstream head(header);
 768         string dummy, label;
 769         head >> dummy;
 770         vector<string> thetaValues(numPartitions, "");
 771         for(int i=0;i<numPartitions;i++){
 772             head >> label >> dummy >> dummy;
 773             thetaValues[i] = label.substr(label.find_last_of('_')+1);
 774         }
 775
 776
 777         vector<double> partitionDiff(numPartitions, 0.0000);
 778
 779         while(referenceFile){
 780
 781             if (m->control_pressed) { break; }
 782
 783             referenceFile >> name >> mean >> lci >> uci;
 784
 785             summaryData tempData;
 786             tempData.name = name;
 787             tempData.refMean = mean;
 788
 789             double difference = 0.0000;
 790
 791             partitionFile >> name;
 792             for(int j=0;j<numPartitions;j++){
 793                 partitionFile >> pMean[j] >> pLCI[j] >> pUCI[j];
 794                 difference += abs(mean - pMean[j]);
 795                 partitionDiff[j] += abs(mean - pMean[j]);;
 796             }
 797
 798             tempData.partMean = pMean;
 799             tempData.partLCI = pLCI;
 800             tempData.partUCI = pUCI;
 801             tempData.difference = difference;
 802             summary.push_back(tempData);
 803
 804             m->gobble(referenceFile);
 805             m->gobble(partitionFile);
 806         }
 807         referenceFile.close();
 808         partitionFile.close();
 809
 810         if (m->control_pressed) { return 0; }
 811
 812         int numOTUs = (int)summary.size();
 813
 814         sort(summary.begin(), summary.end(), summaryFunction);
 815
 816         variables.erase("[tag]");
 817                 string outputFileName = getOutputFileName("parameters", variables);
 818         outputNames.push_back(outputFileName); outputTypes["parameters"].push_back(outputFileName);
 819
 820         ofstream parameterFile;
 821         m->openOutputFile(outputFileName, parameterFile); //((fileRoot + "mix.parameters").c_str());
 822         parameterFile.setf(ios::fixed, ios::floatfield);
 823         parameterFile.setf(ios::showpoint);
 824
 825         double totalDifference =  0.0000;
 826         parameterFile << "Part\tDif2Ref_i\ttheta_i\tpi_i\n";
 827         for(int i=0;i<numPartitions;i++){
 828             if (m->control_pressed) { break; }
 829             parameterFile << i+1 << '\t' << setprecision(2) << partitionDiff[i] << '\t' << thetaValues[i] << '\t' << piValues[i] << endl;
 830             totalDifference += partitionDiff[i];
 831         }
 832         parameterFile.close();
 833
 834         if (m->control_pressed) { return 0; }
 835
 836         string summaryFileName = getOutputFileName("summary", variables);
 837         outputNames.push_back(summaryFileName); outputTypes["summary"].push_back(summaryFileName);
 838
 839         ofstream summaryFile;
 840         m->openOutputFile(summaryFileName, summaryFile); //((fileRoot + "mix.summary").c_str());
 841         summaryFile.setf(ios::fixed, ios::floatfield);
 842         summaryFile.setf(ios::showpoint);
 843
 844
 845         summaryFile << "OTU\tP0.mean";
 846         for(int i=0;i<numPartitions;i++){
 847             summaryFile << "\tP" << i+1 << ".mean\tP" << i+1 << ".lci\tP" << i+1 << ".uci";
 848         }
 849         summaryFile << "\tDifference\tCumFraction" << endl;
 850
 851         double cumDiff = 0.0000;
 852
 853         for(int i=0;i<numOTUs;i++){
 854             if (m->control_pressed) { break; }
 855             summaryFile << summary[i].name << setprecision(2) << '\t' << summary[i].refMean;
 856             for(int j=0;j<numPartitions;j++){
 857                 summaryFile  << '\t' << summary[i].partMean[j] << '\t' << summary[i].partLCI[j] << '\t' << summary[i].partUCI[j];
 858             }
 859
 860             cumDiff += summary[i].difference/totalDifference;
 861             summaryFile << '\t' << summary[i].difference << '\t' << cumDiff << endl;
 862         }
 863         summaryFile.close();
 864
 865         return 0;
 866
 867     }
 868         catch(exception& e) {
 869                 m->errorOut(e, "GetMetaCommunityCommand", "generateSummaryFile");
 870                 exit(1);
 871         }
 872
 873 }
 874 //**********************************************************************************************************************
 875 vector<vector<double> > GetMetaCommunityCommand::generateDistanceMatrix(vector<SharedRAbundVector*>& thisLookup){
 876     try {
 877         vector<vector<double> > results;
 878
 879         Calculator* matrixCalculator;
 880         ValidCalculators validCalculator;
 881         int i = 0;
 882
 883         if (validCalculator.isValidCalculator("matrix", Estimators[i]) == true) {
 884             if (Estimators[i] == "sharedsobs") {
 885                 matrixCalculator = new SharedSobsCS();
 886             }else if (Estimators[i] == "sharedchao") {
 887                 matrixCalculator = new SharedChao1();
 888             }else if (Estimators[i] == "sharedace") {
 889                 matrixCalculator = new SharedAce();
 890             }else if (Estimators[i] == "jabund") {
 891                 matrixCalculator = new JAbund();
 892             }else if (Estimators[i] == "sorabund") {
 893                 matrixCalculator = new SorAbund();
 894             }else if (Estimators[i] == "jclass") {
 895                 matrixCalculator = new Jclass();
 896             }else if (Estimators[i] == "sorclass") {
 897                 matrixCalculator = new SorClass();
 898             }else if (Estimators[i] == "jest") {
 899                 matrixCalculator = new Jest();
 900             }else if (Estimators[i] == "sorest") {
 901                 matrixCalculator = new SorEst();
 902             }else if (Estimators[i] == "thetayc") {
 903                 matrixCalculator = new ThetaYC();
 904             }else if (Estimators[i] == "thetan") {
 905                 matrixCalculator = new ThetaN();
 906             }else if (Estimators[i] == "kstest") {
 907                 matrixCalculator = new KSTest();
 908             }else if (Estimators[i] == "sharednseqs") {
 909                 matrixCalculator = new SharedNSeqs();
 910             }else if (Estimators[i] == "ochiai") {
 911                 matrixCalculator = new Ochiai();
 912             }else if (Estimators[i] == "anderberg") {
 913                 matrixCalculator = new Anderberg();
 914             }else if (Estimators[i] == "kulczynski") {
 915                 matrixCalculator = new Kulczynski();
 916             }else if (Estimators[i] == "kulczynskicody") {
 917                 matrixCalculator = new KulczynskiCody();
 918             }else if (Estimators[i] == "lennon") {
 919                 matrixCalculator = new Lennon();
 920             }else if (Estimators[i] == "morisitahorn") {
 921                 matrixCalculator = new MorHorn();
 922             }else if (Estimators[i] == "braycurtis") {
 923                 matrixCalculator = new BrayCurtis();
 924             }else if (Estimators[i] == "whittaker") {
 925                 matrixCalculator = new Whittaker();
 926             }else if (Estimators[i] == "odum") {
 927                 matrixCalculator = new Odum();
 928             }else if (Estimators[i] == "canberra") {
 929                 matrixCalculator = new Canberra();
 930             }else if (Estimators[i] == "structeuclidean") {
 931                 matrixCalculator = new StructEuclidean();
 932             }else if (Estimators[i] == "structchord") {
 933                 matrixCalculator = new StructChord();
 934             }else if (Estimators[i] == "hellinger") {
 935                 matrixCalculator = new Hellinger();
 936             }else if (Estimators[i] == "manhattan") {
 937                 matrixCalculator = new Manhattan();
 938             }else if (Estimators[i] == "structpearson") {
 939                 matrixCalculator = new StructPearson();
 940             }else if (Estimators[i] == "soergel") {
 941                 matrixCalculator = new Soergel();
 942             }else if (Estimators[i] == "spearman") {
 943                 matrixCalculator = new Spearman();
 944             }else if (Estimators[i] == "structkulczynski") {
 945                 matrixCalculator = new StructKulczynski();
 946             }else if (Estimators[i] == "speciesprofile") {
 947                 matrixCalculator = new SpeciesProfile();
 948             }else if (Estimators[i] == "hamming") {
 949                 matrixCalculator = new Hamming();
 950             }else if (Estimators[i] == "structchi2") {
 951                 matrixCalculator = new StructChi2();
 952             }else if (Estimators[i] == "gower") {
 953                 matrixCalculator = new Gower();
 954             }else if (Estimators[i] == "memchi2") {
 955                 matrixCalculator = new MemChi2();
 956             }else if (Estimators[i] == "memchord") {
 957                 matrixCalculator = new MemChord();
 958             }else if (Estimators[i] == "memeuclidean") {
 959                 matrixCalculator = new MemEuclidean();
 960             }else if (Estimators[i] == "mempearson") {
 961                 matrixCalculator = new MemPearson();
 962             }else if (Estimators[i] == "jsd") {
 963                 matrixCalculator = new JSD();
 964             }else if (Estimators[i] == "rjsd") {
 965                 matrixCalculator = new RJSD();
 966             }else {
 967                 m->mothurOut("[ERROR]: " + Estimators[i] + " is not a valid calculator, please correct.\n"); m->control_pressed = true; return results;
 968             }
 969         }
 970
 971         //calc distances
 972         vector< vector< vector<seqDist> > > calcDistsTotals;  //each iter, then each groupCombos dists. this will be used to make .dist files
 973         vector< vector<seqDist> > calcDists; calcDists.resize(1);
 974
 975         for (int thisIter = 0; thisIter < iters+1; thisIter++) {
 976
 977             vector<SharedRAbundVector*> thisItersLookup = thisLookup;
 978
 979             if (subsample && (thisIter != 0)) {
 980                 SubSample sample;
 981                 vector<string> tempLabels; //dont need since we arent printing the sampled sharedRabunds
 982
 983                 //make copy of lookup so we don't get access violations
 984                 vector<SharedRAbundVector*> newLookup;
 985                 for (int k = 0; k < thisItersLookup.size(); k++) {
 986                     SharedRAbundVector* temp = new SharedRAbundVector();
 987                     temp->setLabel(thisItersLookup[k]->getLabel());
 988                     temp->setGroup(thisItersLookup[k]->getGroup());
 989                     newLookup.push_back(temp);
 990                 }
 991
 992                 //for each bin
 993                 for (int k = 0; k < thisItersLookup[0]->getNumBins(); k++) {
 994                     if (m->control_pressed) { for (int j = 0; j < newLookup.size(); j++) {  delete newLookup[j];  } return results; }
 995                     for (int j = 0; j < thisItersLookup.size(); j++) { newLookup[j]->push_back(thisItersLookup[j]->getAbundance(k), thisItersLookup[j]->getGroup()); }
 996                 }
 997
 998                 tempLabels = sample.getSample(newLookup, subsampleSize);
 999                 thisItersLookup = newLookup;
1000             }
1001
1002
1003             driver(thisItersLookup, calcDists, matrixCalculator);
1004
1005             if (subsample && (thisIter != 0)) {
1006                 if((thisIter) % 100 == 0){      m->mothurOutJustToScreen(toString(thisIter)+"\n");              }
1007                 calcDistsTotals.push_back(calcDists);
1008                 for (int i = 0; i < calcDists.size(); i++) {
1009                     for (int j = 0; j < calcDists[i].size(); j++) {
1010                         if (m->debug) {  m->mothurOut("[DEBUG]: Results: iter = " + toString(thisIter) + ", " + thisLookup[calcDists[i][j].seq1]->getGroup() + " - " + thisLookup[calcDists[i][j].seq2]->getGroup() + " distance = " + toString(calcDists[i][j].dist) + ".\n");  }
1011                     }
1012                 }
1013                 //clean up memory
1014                 for (int i = 0; i < thisItersLookup.size(); i++) { delete thisItersLookup[i]; }
1015                 thisItersLookup.clear();
1016             }else { //print results for whole dataset
1017                 for (int i = 0; i < calcDists.size(); i++) {
1018                     if (m->control_pressed) { break; }
1019
1020                     //initialize matrix
1021                     results.resize(thisLookup.size());
1022                     for (int k = 0; k < thisLookup.size(); k++) {  results[k].resize(thisLookup.size(), 0.0); }
1023
1024                     for (int j = 0; j < calcDists[i].size(); j++) {
1025                         int row = calcDists[i][j].seq1;
1026                         int column = calcDists[i][j].seq2;
1027                         double dist = calcDists[i][j].dist;
1028
1029                         results[row][column] = dist;
1030                         results[column][row] = dist;
1031                     }
1032                 }
1033             }
1034             for (int i = 0; i < calcDists.size(); i++) {  calcDists[i].clear(); }
1035                 }
1036
1037         if (iters != 0) {
1038             //we need to find the average distance and standard deviation for each groups distance
1039             vector< vector<seqDist>  > calcAverages = m->getAverages(calcDistsTotals, "average");
1040
1041             //print results
1042             for (int i = 0; i < calcDists.size(); i++) {
1043                 results.resize(thisLookup.size());
1044                 for (int k = 0; k < thisLookup.size(); k++) {  results[k].resize(thisLookup.size(), 0.0); }
1045
1046                 for (int j = 0; j < calcAverages[i].size(); j++) {
1047                     int row = calcAverages[i][j].seq1;
1048                     int column = calcAverages[i][j].seq2;
1049                     float dist = calcAverages[i][j].dist;
1050
1051                     results[row][column] = dist;
1052                     results[column][row] = dist;
1053                 }
1054             }
1055         }
1056
1057
1058         return results;
1059     }
1060     catch(exception& e) {
1061         m->errorOut(e, "GetMetaCommunityCommand", "generateDistanceMatrix");
1062         exit(1);
1063     }
1064 }
1065 /**************************************************************************************************/
1066 int GetMetaCommunityCommand::driver(vector<SharedRAbundVector*> thisLookup, vector< vector<seqDist> >& calcDists, Calculator* matrixCalculator) {
1067         try {
1068                 vector<SharedRAbundVector*> subset;
1069
1070                 for (int k = 0; k < thisLookup.size(); k++) { // pass cdd each set of groups to compare
1071
1072                         for (int l = 0; l < k; l++) {
1073
1074                                 if (k != l) { //we dont need to similiarity of a groups to itself
1075                                         subset.clear(); //clear out old pair of sharedrabunds
1076                                         //add new pair of sharedrabunds
1077                                         subset.push_back(thisLookup[k]); subset.push_back(thisLookup[l]);
1078
1079
1080
1081                     //if this calc needs all groups to calculate the pair load all groups
1082                     if (matrixCalculator->getNeedsAll()) {
1083                         //load subset with rest of lookup for those calcs that need everyone to calc for a pair
1084                         for (int w = 0; w < thisLookup.size(); w++) {
1085                             if ((w != k) && (w != l)) { subset.push_back(thisLookup[w]); }
1086                         }
1087                     }
1088
1089                     vector<double> tempdata = matrixCalculator->getValues(subset); //saves the calculator outputs
1090
1091                     if (m->control_pressed) { return 1; }
1092
1093                     seqDist temp(l, k, tempdata[0]);
1094                     //cout << l << '\t' << k << '\t' <<  tempdata[0] << endl;
1095                     calcDists[0].push_back(temp);
1096                 }
1097
1098                         }
1099                 }
1100
1101                 return 0;
1102         }
1103         catch(exception& e) {
1104                 m->errorOut(e, "MatrixOutputCommand", "driver");
1105                 exit(1);
1106         }
1107 }
1108 //**********************************************************************************************************************
1109