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],[method],mix.fit"; }
  70         else if (type == "relabund")    {  pattern = "[filename],[distance],[method],[tag],mix.relabund"; }
  71         else if (type == "design")      {  pattern = "[filename],[distance],[method],mix.design"; }
  72         else if (type == "matrix")      {  pattern = "[filename],[distance],[method],[tag],mix.posterior"; }
  73         else if (type == "parameters")  {  pattern = "[filename],[distance],[method],mix.parameters"; }
  74         else if (type == "summary")  {  pattern = "[filename],[distance],[method],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         variables["[method]"] = method;
 370                 string outputFileName = getOutputFileName("fit", variables);
 371         outputNames.push_back(outputFileName); outputTypes["fit"].push_back(outputFileName);
 372
 373                 //divide the partitions between the processors
 374                 vector< vector<int> > dividedPartitions;
 375         vector< vector<string> > rels, matrix;
 376         vector<string> doneFlags;
 377         dividedPartitions.resize(processors);
 378         rels.resize(processors);
 379         matrix.resize(processors);
 380
 381         //for each file group figure out which process will complete it
 382         //want to divide the load intelligently so the big files are spread between processes
 383         for (int i=1; i<=maxpartitions; i++) {
 384             //cout << i << endl;
 385             int processToAssign = (i+1) % processors;
 386             if (processToAssign == 0) { processToAssign = processors; }
 387
 388             if (m->debug) { m->mothurOut("[DEBUG]: assigning " + toString(i) + " to process " + toString(processToAssign-1) + "\n"); }
 389             dividedPartitions[(processToAssign-1)].push_back(i);
 390
 391             variables["[tag]"] = toString(i);
 392             string relName = getOutputFileName("relabund", variables);
 393             string mName = getOutputFileName("matrix", variables);
 394             rels[(processToAssign-1)].push_back(relName);
 395             matrix[(processToAssign-1)].push_back(mName);
 396         }
 397
 398         for (int i = 0; i < processors; i++) { //read from everyone elses, just write to yours
 399             string tempDoneFile = m->getRootName(m->getSimpleName(sharedfile)) + toString(i) + ".done.temp";
 400             doneFlags.push_back(tempDoneFile);
 401             ofstream out;
 402             m->openOutputFile(tempDoneFile, out); //clear out
 403             out.close();
 404         }
 405
 406
 407 #if defined (__APPLE__) || (__MACH__) || (linux) || (__linux) || (__linux__) || (__unix__) || (__unix)
 408
 409                 //loop through and create all the processes you want
 410                 while (process != processors) {
 411                         int pid = fork();
 412
 413                         if (pid > 0) {
 414                                 processIDS.push_back(pid);  //create map from line number to pid so you can append files in correct order later
 415                                 process++;
 416                         }else if (pid == 0){
 417                 outputNames.clear();
 418                                 num = processDriver(thislookup, dividedPartitions[process], (outputFileName + toString(getpid())), rels[process], matrix[process], doneFlags, process);
 419
 420                 //pass numSeqs to parent
 421                                 ofstream out;
 422                                 string tempFile = toString(getpid()) + ".outputNames.temp";
 423                                 m->openOutputFile(tempFile, out);
 424                 out << num << endl;
 425                 out << outputNames.size() << endl;
 426                                 for (int i = 0; i < outputNames.size(); i++) { out << outputNames[i] << endl; }
 427                                 out.close();
 428
 429                                 exit(0);
 430                         }else {
 431                                 m->mothurOut("[ERROR]: unable to spawn the necessary processes."); m->mothurOutEndLine();
 432                                 for (int i = 0; i < processIDS.size(); i++) { kill (processIDS[i], SIGINT); }
 433                                 exit(0);
 434                         }
 435                 }
 436
 437                 //do my part
 438         if (method == "dmm") {  m->mothurOut("K\tNLE\t\tlogDet\tBIC\t\tAIC\t\tLaplace\n");  }
 439         else {
 440             m->mothurOut("K\tCH\t");
 441             for (int i = 0; i < thislookup.size(); i++) {  m->mothurOut(thislookup[i]->getGroup() + '\t'); }
 442             m->mothurOut("\n");
 443         }
 444                 minPartition = processDriver(thislookup, dividedPartitions[0], outputFileName, rels[0], matrix[0], doneFlags, 0);
 445
 446                 //force parent to wait until all the processes are done
 447                 for (int i=0;i<processIDS.size();i++) {
 448                         int temp = processIDS[i];
 449                         wait(&temp);
 450                 }
 451
 452         vector<string> tempOutputNames = outputNames;
 453         for (int i=0;i<processIDS.size();i++) {
 454             ifstream in;
 455                         string tempFile = toString(processIDS[i]) + ".outputNames.temp";
 456                         m->openInputFile(tempFile, in);
 457                         if (!in.eof()) {
 458                 int tempNum = 0;
 459                 in >> tempNum; m->gobble(in);
 460                 if (tempNum < minPartition) { minPartition = tempNum; }
 461                 in >> tempNum; m->gobble(in);
 462                 for (int i = 0; i < tempNum; i++) {
 463                     string tempName = "";
 464                     in >> tempName; m->gobble(in);
 465                     tempOutputNames.push_back(tempName);
 466                 }
 467             }
 468                         in.close(); m->mothurRemove(tempFile);
 469
 470             m->appendFilesWithoutHeaders(outputFileName + toString(processIDS[i]), outputFileName);
 471             m->mothurRemove(outputFileName + toString(processIDS[i]));
 472         }
 473
 474         if (processors > 1) {
 475             outputNames.clear();
 476             for (int i = 0; i < tempOutputNames.size(); i++) { //remove files if needed
 477                 string name = tempOutputNames[i];
 478                 vector<string> parts;
 479                 m->splitAtChar(name, parts, '.');
 480                 bool keep = true;
 481                 if (((parts[parts.size()-1] == "relabund") || (parts[parts.size()-1] == "posterior")) && (parts[parts.size()-2] == "mix")) {
 482                     string tempNum = parts[parts.size()-3];
 483                     int num;  m->mothurConvert(tempNum, num);
 484                     //if (num > minPartition) {
 485                      //   m->mothurRemove(tempOutputNames[i]);
 486                     //    keep = false; if (m->debug) { m->mothurOut("[DEBUG]: removing " + tempOutputNames[i] + ".\n"); }
 487                     //}
 488                 }
 489                 if (keep) { outputNames.push_back(tempOutputNames[i]); }
 490             }
 491
 492             //reorder fit file
 493             ifstream in;
 494             m->openInputFile(outputFileName, in);
 495             string headers = m->getline(in); m->gobble(in);
 496
 497             map<int, string> file;
 498             while (!in.eof()) {
 499                 string numString, line;
 500                 int num;
 501                 in >> numString; line = m->getline(in); m->gobble(in);
 502                 m->mothurConvert(numString, num);
 503                 file[num] = line;
 504             }
 505             in.close();
 506             ofstream out;
 507             m->openOutputFile(outputFileName, out);
 508             out << headers << endl;
 509             for (map<int, string>::iterator it = file.begin(); it != file.end(); it++) {
 510                 out << it->first << '\t' << it->second << endl;
 511                 if (m->debug) { m->mothurOut("[DEBUG]: printing: " + toString(it->first) + '\t' + it->second + ".\n"); }
 512             }
 513             out.close();
 514         }
 515
 516 #else
 517         m->mothurOut("K\tNLE\t\tlogDet\tBIC\t\tAIC\t\tLaplace\n");
 518                 minPartition = processDriver(thislookup, dividedPartitions[0], outputFileName, rels[0], matrix[0], doneFlags, 0);
 519 #endif
 520         for (int i = 0; i < processors; i++) { //read from everyone elses, just write to yours
 521             string tempDoneFile = m->getRootName(m->getSimpleName(sharedfile)) + toString(i) + ".done.temp";
 522             m->mothurRemove(tempDoneFile);
 523         }
 524
 525         if (m->control_pressed) { return 0; }
 526
 527         if (m->debug) { m->mothurOut("[DEBUG]: minPartition = " + toString(minPartition) + "\n"); }
 528
 529         //run generate Summary function for smallest minPartition
 530         variables["[tag]"] = toString(minPartition);
 531         vector<double> piValues = generateDesignFile(minPartition, variables);
 532         if (method == "dmm") {  generateSummaryFile(minPartition, variables, piValues); } //pam doesn't make a relabund file
 533
 534         return 0;
 535
 536     }
 537         catch(exception& e) {
 538                 m->errorOut(e, "GetMetaCommunityCommand", "createProcesses");
 539                 exit(1);
 540         }
 541 }
 542 //**********************************************************************************************************************
 543 int GetMetaCommunityCommand::processDriver(vector<SharedRAbundVector*>& thislookup, vector<int>& parts, string outputFileName, vector<string> relabunds, vector<string> matrix, vector<string> doneFlags, int processID){
 544         try {
 545
 546         double minLaplace = 1e10;
 547         int minPartition = 1;
 548         vector<double> minSilhouettes; minSilhouettes.resize(thislookup.size(), 0);
 549
 550                 ofstream fitData, silData;
 551         if (method == "dmm") {
 552             m->openOutputFile(outputFileName, fitData);
 553             fitData.setf(ios::fixed, ios::floatfield);
 554             fitData.setf(ios::showpoint);
 555             fitData << "K\tNLE\tlogDet\tBIC\tAIC\tLaplace" << endl;
 556         }else if((method == "pam") || (method == "kmeans")) { //because ch is looking of maximal value
 557             minLaplace = 0;
 558             m->openOutputFile(outputFileName, silData);
 559             silData.setf(ios::fixed, ios::floatfield);
 560             silData.setf(ios::showpoint);
 561             silData << "K\tCH\t";
 562             for (int i = 0; i < thislookup.size(); i++) { silData << thislookup[i]->getGroup() << '\t';  }
 563             silData << endl;
 564         }
 565
 566         cout.setf(ios::fixed, ios::floatfield);
 567         cout.setf(ios::showpoint);
 568
 569         vector< vector<int> > sharedMatrix;
 570         vector<string> thisGroups;
 571         for (int i = 0; i < thislookup.size(); i++) { sharedMatrix.push_back(thislookup[i]->getAbundances()); thisGroups.push_back(thislookup[i]->getGroup()); }
 572
 573         vector< vector<double> > dists; //do we want to output this matrix??
 574         if ((method == "pam") || (method == "kmeans")) {  dists = generateDistanceMatrix(thislookup);  }
 575
 576         if (m->debug) {
 577             m->mothurOut("[DEBUG]: dists = \n");
 578             for (int i = 0; i < dists.size(); i++) {
 579                 if (m->control_pressed) { break; }
 580                 m->mothurOut("[DEBUG]: i = " + toString(i) + '\t');
 581                 for (int j = 0; j < i; j++) { m->mothurOut(toString(dists[i][j]) +"\t"); }
 582                 m->mothurOut("\n");
 583             }
 584         }
 585
 586         for(int i=0;i<parts.size();i++){
 587
 588             int numPartitions = parts[i];
 589
 590             if (m->debug) { m->mothurOut("[DEBUG]: running partition " + toString(numPartitions) + "\n"); }
 591
 592             if (m->control_pressed) { break; }
 593
 594             //check to see if anyone else is done
 595             for (int j = 0; j < doneFlags.size(); j++) {
 596                 if (!m->isBlank(doneFlags[j])) { //another process has finished
 597                     //are they done at a lower partition?
 598                     ifstream in;
 599                     m->openInputFile(doneFlags[j], in);
 600                     int tempNum;
 601                     in >> tempNum; in.close();
 602                     if (tempNum < numPartitions) { break; } //quit, because someone else has finished
 603                 }
 604             }
 605
 606             CommunityTypeFinder* finder = NULL;
 607             if (method == "dmm")            {   finder = new qFinderDMM(sharedMatrix, numPartitions);   }
 608             else if (method == "kmeans")    {   finder = new KMeans(sharedMatrix, numPartitions);       }
 609             else if (method == "pam")       {   finder = new Pam(sharedMatrix, dists, numPartitions);                 }
 610             else {
 611                 if (i == 0) {  m->mothurOut(method + " is not a valid method option. I will run the command using dmm.\n"); }
 612                 finder = new qFinderDMM(sharedMatrix, numPartitions);
 613             }
 614
 615             string relabund = relabunds[i];
 616             string matrixName = matrix[i];
 617             outputNames.push_back(matrixName); outputTypes["matrix"].push_back(matrixName);
 618
 619             finder->printZMatrix(matrixName, thisGroups);
 620
 621             double chi; vector<double> silhouettes;
 622             if (method == "dmm") {
 623                 double laplace = finder->getLaplace();
 624                 if(laplace < minLaplace){
 625                     minPartition = numPartitions;
 626                     minLaplace = laplace;
 627                 }
 628             }else {
 629                 chi = finder->calcCHIndex(dists);
 630                 silhouettes = finder->calcSilhouettes(dists);
 631                 if (chi > minLaplace) { //save partition with maximum ch index score
 632                     minPartition = numPartitions;
 633                     minLaplace = chi;
 634                     minSilhouettes = silhouettes;
 635                 }
 636             }
 637
 638             if (method == "dmm") {
 639                 finder->printFitData(cout, minLaplace);
 640                 finder->printFitData(fitData);
 641                 finder->printRelAbund(relabund, m->currentSharedBinLabels);
 642                 outputNames.push_back(relabund); outputTypes["relabund"].push_back(relabund);
 643             }else if ((method == "pam") || (method == "kmeans")) { //print silouettes and ch values
 644                 finder->printSilData(cout, chi, silhouettes);
 645                 finder->printSilData(silData, chi, silhouettes);
 646                 if (method == "kmeans") {
 647                     finder->printRelAbund(relabund, m->currentSharedBinLabels);
 648                     outputNames.push_back(relabund); outputTypes["relabund"].push_back(relabund);
 649                 }
 650             }
 651             delete finder;
 652
 653             if(optimizegap != -1 && (numPartitions - minPartition) >= optimizegap && numPartitions >= minpartitions){
 654                 string tempDoneFile = m->getRootName(m->getSimpleName(sharedfile)) + toString(processID) + ".done.temp";
 655                 ofstream outDone;
 656                 m->openOutputFile(tempDoneFile, outDone);
 657                 outDone << minPartition << endl;
 658                 outDone.close();
 659                 break;
 660             }
 661         }
 662         if (method == "dmm") { fitData.close(); }
 663
 664         if (m->control_pressed) { return 0; }
 665
 666         return minPartition;
 667     }
 668         catch(exception& e) {
 669                 m->errorOut(e, "GetMetaCommunityCommand", "processDriver");
 670                 exit(1);
 671         }
 672 }
 673 /**************************************************************************************************/
 674
 675 vector<double> GetMetaCommunityCommand::generateDesignFile(int numPartitions, map<string,string> variables){
 676     try {
 677         vector<double> piValues(numPartitions, 0);
 678
 679         ifstream postFile;
 680         variables["[tag]"] = toString(numPartitions);
 681         string input = getOutputFileName("matrix", variables);
 682         m->openInputFile(input, postFile);//((fileRoot + toString(numPartitions) + "mix.posterior").c_str()); //matrix file
 683         variables.erase("[tag]");
 684                 string outputFileName = getOutputFileName("design", variables);
 685         ofstream designFile;
 686         m->openOutputFile(outputFileName, designFile);
 687         outputNames.push_back(outputFileName); outputTypes["design"].push_back(outputFileName);
 688
 689
 690         vector<string> titles(numPartitions);
 691
 692         for(int i=0;i<numPartitions;i++){   postFile >> titles[i];  }
 693
 694         double posterior;
 695         string sampleName;
 696         int numSamples = 0;
 697
 698         while(postFile){
 699
 700             if (m->control_pressed) { break; }
 701
 702             double maxPosterior = 0.0000;
 703             int maxPartition = -1;
 704
 705             postFile >> sampleName;
 706
 707             for(int i=0;i<numPartitions;i++){
 708
 709                 postFile >> posterior;
 710                 if(posterior > maxPosterior){
 711                     maxPosterior = posterior;
 712                     maxPartition = i;
 713                 }
 714                 piValues[i] += posterior;
 715
 716             }
 717
 718             designFile << sampleName << '\t' << titles[maxPartition] << endl;
 719
 720             numSamples++;
 721             m->gobble(postFile);
 722         }
 723         for(int i=0;i<numPartitions;i++){
 724             piValues[i] /= (double)numSamples;
 725         }
 726
 727
 728         postFile.close();
 729         designFile.close();
 730
 731         return piValues;
 732     }
 733         catch(exception& e) {
 734                 m->errorOut(e, "GetMetaCommunityCommand", "generateDesignFile");
 735                 exit(1);
 736         }
 737 }
 738
 739 /**************************************************************************************************/
 740
 741 inline bool summaryFunction(summaryData i, summaryData j){ return i.difference > j.difference;   }
 742
 743 /**************************************************************************************************/
 744 int GetMetaCommunityCommand::generateSummaryFile(int numPartitions, map<string,string> v, vector<double> piValues){
 745     try {
 746         vector<summaryData> summary;
 747
 748         vector<double> pMean(numPartitions, 0);
 749         vector<double> pLCI(numPartitions, 0);
 750         vector<double> pUCI(numPartitions, 0);
 751
 752         string name, header;
 753         double mean, lci, uci;
 754
 755         ifstream referenceFile;
 756         map<string, string> variables;
 757         variables["[filename]"] = v["[filename]"];
 758         variables["[distance]"] = v["[distance]"];
 759         variables["]method]"] = method;
 760         variables["[tag]"] = "1";
 761         string reference = getOutputFileName("relabund", variables);
 762         m->openInputFile(reference, referenceFile); //((fileRoot + label + ".1mix.relabund").c_str());
 763         variables["[tag]"] = toString(numPartitions);
 764         string partFile = getOutputFileName("relabund", variables);
 765         ifstream partitionFile;
 766         m->openInputFile(partFile, partitionFile); //((fileRoot + toString(numPartitions) + "mix.relabund").c_str());
 767
 768         header = m->getline(referenceFile);
 769         header = m->getline(partitionFile);
 770         stringstream head(header);
 771         string dummy, label;
 772         head >> dummy;
 773         vector<string> thetaValues(numPartitions, "");
 774         for(int i=0;i<numPartitions;i++){
 775             head >> label >> dummy >> dummy;
 776             thetaValues[i] = label.substr(label.find_last_of('_')+1);
 777         }
 778
 779
 780         vector<double> partitionDiff(numPartitions, 0.0000);
 781
 782         while(referenceFile){
 783
 784             if (m->control_pressed) { break; }
 785
 786             referenceFile >> name >> mean >> lci >> uci;
 787
 788             summaryData tempData;
 789             tempData.name = name;
 790             tempData.refMean = mean;
 791
 792             double difference = 0.0000;
 793
 794             partitionFile >> name;
 795             for(int j=0;j<numPartitions;j++){
 796                 partitionFile >> pMean[j] >> pLCI[j] >> pUCI[j];
 797                 difference += abs(mean - pMean[j]);
 798                 partitionDiff[j] += abs(mean - pMean[j]);;
 799             }
 800
 801             tempData.partMean = pMean;
 802             tempData.partLCI = pLCI;
 803             tempData.partUCI = pUCI;
 804             tempData.difference = difference;
 805             summary.push_back(tempData);
 806
 807             m->gobble(referenceFile);
 808             m->gobble(partitionFile);
 809         }
 810         referenceFile.close();
 811         partitionFile.close();
 812
 813         if (m->control_pressed) { return 0; }
 814
 815         int numOTUs = (int)summary.size();
 816
 817         sort(summary.begin(), summary.end(), summaryFunction);
 818
 819         variables.erase("[tag]");
 820                 string outputFileName = getOutputFileName("parameters", variables);
 821         outputNames.push_back(outputFileName); outputTypes["parameters"].push_back(outputFileName);
 822
 823         ofstream parameterFile;
 824         m->openOutputFile(outputFileName, parameterFile); //((fileRoot + "mix.parameters").c_str());
 825         parameterFile.setf(ios::fixed, ios::floatfield);
 826         parameterFile.setf(ios::showpoint);
 827
 828         double totalDifference =  0.0000;
 829         parameterFile << "Part\tDif2Ref_i\ttheta_i\tpi_i\n";
 830         for(int i=0;i<numPartitions;i++){
 831             if (m->control_pressed) { break; }
 832             parameterFile << i+1 << '\t' << setprecision(2) << partitionDiff[i] << '\t' << thetaValues[i] << '\t' << piValues[i] << endl;
 833             totalDifference += partitionDiff[i];
 834         }
 835         parameterFile.close();
 836
 837         if (m->control_pressed) { return 0; }
 838
 839         string summaryFileName = getOutputFileName("summary", variables);
 840         outputNames.push_back(summaryFileName); outputTypes["summary"].push_back(summaryFileName);
 841
 842         ofstream summaryFile;
 843         m->openOutputFile(summaryFileName, summaryFile); //((fileRoot + "mix.summary").c_str());
 844         summaryFile.setf(ios::fixed, ios::floatfield);
 845         summaryFile.setf(ios::showpoint);
 846
 847
 848         summaryFile << "OTU\tP0.mean";
 849         for(int i=0;i<numPartitions;i++){
 850             summaryFile << "\tP" << i+1 << ".mean\tP" << i+1 << ".lci\tP" << i+1 << ".uci";
 851         }
 852         summaryFile << "\tDifference\tCumFraction" << endl;
 853
 854         double cumDiff = 0.0000;
 855
 856         for(int i=0;i<numOTUs;i++){
 857             if (m->control_pressed) { break; }
 858             summaryFile << summary[i].name << setprecision(2) << '\t' << summary[i].refMean;
 859             for(int j=0;j<numPartitions;j++){
 860                 summaryFile  << '\t' << summary[i].partMean[j] << '\t' << summary[i].partLCI[j] << '\t' << summary[i].partUCI[j];
 861             }
 862
 863             cumDiff += summary[i].difference/totalDifference;
 864             summaryFile << '\t' << summary[i].difference << '\t' << cumDiff << endl;
 865         }
 866         summaryFile.close();
 867
 868         return 0;
 869
 870     }
 871         catch(exception& e) {
 872                 m->errorOut(e, "GetMetaCommunityCommand", "generateSummaryFile");
 873                 exit(1);
 874         }
 875
 876 }
 877 //**********************************************************************************************************************
 878 vector<vector<double> > GetMetaCommunityCommand::generateDistanceMatrix(vector<SharedRAbundVector*>& thisLookup){
 879     try {
 880         vector<vector<double> > results;
 881
 882         Calculator* matrixCalculator;
 883         ValidCalculators validCalculator;
 884         int i = 0;
 885
 886         if (validCalculator.isValidCalculator("matrix", Estimators[i]) == true) {
 887             if (Estimators[i] == "sharedsobs") {
 888                 matrixCalculator = new SharedSobsCS();
 889             }else if (Estimators[i] == "sharedchao") {
 890                 matrixCalculator = new SharedChao1();
 891             }else if (Estimators[i] == "sharedace") {
 892                 matrixCalculator = new SharedAce();
 893             }else if (Estimators[i] == "jabund") {
 894                 matrixCalculator = new JAbund();
 895             }else if (Estimators[i] == "sorabund") {
 896                 matrixCalculator = new SorAbund();
 897             }else if (Estimators[i] == "jclass") {
 898                 matrixCalculator = new Jclass();
 899             }else if (Estimators[i] == "sorclass") {
 900                 matrixCalculator = new SorClass();
 901             }else if (Estimators[i] == "jest") {
 902                 matrixCalculator = new Jest();
 903             }else if (Estimators[i] == "sorest") {
 904                 matrixCalculator = new SorEst();
 905             }else if (Estimators[i] == "thetayc") {
 906                 matrixCalculator = new ThetaYC();
 907             }else if (Estimators[i] == "thetan") {
 908                 matrixCalculator = new ThetaN();
 909             }else if (Estimators[i] == "kstest") {
 910                 matrixCalculator = new KSTest();
 911             }else if (Estimators[i] == "sharednseqs") {
 912                 matrixCalculator = new SharedNSeqs();
 913             }else if (Estimators[i] == "ochiai") {
 914                 matrixCalculator = new Ochiai();
 915             }else if (Estimators[i] == "anderberg") {
 916                 matrixCalculator = new Anderberg();
 917             }else if (Estimators[i] == "kulczynski") {
 918                 matrixCalculator = new Kulczynski();
 919             }else if (Estimators[i] == "kulczynskicody") {
 920                 matrixCalculator = new KulczynskiCody();
 921             }else if (Estimators[i] == "lennon") {
 922                 matrixCalculator = new Lennon();
 923             }else if (Estimators[i] == "morisitahorn") {
 924                 matrixCalculator = new MorHorn();
 925             }else if (Estimators[i] == "braycurtis") {
 926                 matrixCalculator = new BrayCurtis();
 927             }else if (Estimators[i] == "whittaker") {
 928                 matrixCalculator = new Whittaker();
 929             }else if (Estimators[i] == "odum") {
 930                 matrixCalculator = new Odum();
 931             }else if (Estimators[i] == "canberra") {
 932                 matrixCalculator = new Canberra();
 933             }else if (Estimators[i] == "structeuclidean") {
 934                 matrixCalculator = new StructEuclidean();
 935             }else if (Estimators[i] == "structchord") {
 936                 matrixCalculator = new StructChord();
 937             }else if (Estimators[i] == "hellinger") {
 938                 matrixCalculator = new Hellinger();
 939             }else if (Estimators[i] == "manhattan") {
 940                 matrixCalculator = new Manhattan();
 941             }else if (Estimators[i] == "structpearson") {
 942                 matrixCalculator = new StructPearson();
 943             }else if (Estimators[i] == "soergel") {
 944                 matrixCalculator = new Soergel();
 945             }else if (Estimators[i] == "spearman") {
 946                 matrixCalculator = new Spearman();
 947             }else if (Estimators[i] == "structkulczynski") {
 948                 matrixCalculator = new StructKulczynski();
 949             }else if (Estimators[i] == "speciesprofile") {
 950                 matrixCalculator = new SpeciesProfile();
 951             }else if (Estimators[i] == "hamming") {
 952                 matrixCalculator = new Hamming();
 953             }else if (Estimators[i] == "structchi2") {
 954                 matrixCalculator = new StructChi2();
 955             }else if (Estimators[i] == "gower") {
 956                 matrixCalculator = new Gower();
 957             }else if (Estimators[i] == "memchi2") {
 958                 matrixCalculator = new MemChi2();
 959             }else if (Estimators[i] == "memchord") {
 960                 matrixCalculator = new MemChord();
 961             }else if (Estimators[i] == "memeuclidean") {
 962                 matrixCalculator = new MemEuclidean();
 963             }else if (Estimators[i] == "mempearson") {
 964                 matrixCalculator = new MemPearson();
 965             }else if (Estimators[i] == "jsd") {
 966                 matrixCalculator = new JSD();
 967             }else if (Estimators[i] == "rjsd") {
 968                 matrixCalculator = new RJSD();
 969             }else {
 970                 m->mothurOut("[ERROR]: " + Estimators[i] + " is not a valid calculator, please correct.\n"); m->control_pressed = true; return results;
 971             }
 972         }
 973
 974         //calc distances
 975         vector< vector< vector<seqDist> > > calcDistsTotals;  //each iter, then each groupCombos dists. this will be used to make .dist files
 976         vector< vector<seqDist> > calcDists; calcDists.resize(1);
 977
 978         for (int thisIter = 0; thisIter < iters+1; thisIter++) {
 979
 980             vector<SharedRAbundVector*> thisItersLookup = thisLookup;
 981
 982             if (subsample && (thisIter != 0)) {
 983                 SubSample sample;
 984                 vector<string> tempLabels; //dont need since we arent printing the sampled sharedRabunds
 985
 986                 //make copy of lookup so we don't get access violations
 987                 vector<SharedRAbundVector*> newLookup;
 988                 for (int k = 0; k < thisItersLookup.size(); k++) {
 989                     SharedRAbundVector* temp = new SharedRAbundVector();
 990                     temp->setLabel(thisItersLookup[k]->getLabel());
 991                     temp->setGroup(thisItersLookup[k]->getGroup());
 992                     newLookup.push_back(temp);
 993                 }
 994
 995                 //for each bin
 996                 for (int k = 0; k < thisItersLookup[0]->getNumBins(); k++) {
 997                     if (m->control_pressed) { for (int j = 0; j < newLookup.size(); j++) {  delete newLookup[j];  } return results; }
 998                     for (int j = 0; j < thisItersLookup.size(); j++) { newLookup[j]->push_back(thisItersLookup[j]->getAbundance(k), thisItersLookup[j]->getGroup()); }
 999                 }
1000
1001                 tempLabels = sample.getSample(newLookup, subsampleSize);
1002                 thisItersLookup = newLookup;
1003             }
1004
1005
1006             driver(thisItersLookup, calcDists, matrixCalculator);
1007
1008             if (subsample && (thisIter != 0)) {
1009                 if((thisIter) % 100 == 0){      m->mothurOutJustToScreen(toString(thisIter)+"\n");              }
1010                 calcDistsTotals.push_back(calcDists);
1011                 for (int i = 0; i < calcDists.size(); i++) {
1012                     for (int j = 0; j < calcDists[i].size(); j++) {
1013                         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");  }
1014                     }
1015                 }
1016                 //clean up memory
1017                 for (int i = 0; i < thisItersLookup.size(); i++) { delete thisItersLookup[i]; }
1018                 thisItersLookup.clear();
1019             }else { //print results for whole dataset
1020                 for (int i = 0; i < calcDists.size(); i++) {
1021                     if (m->control_pressed) { break; }
1022
1023                     //initialize matrix
1024                     results.resize(thisLookup.size());
1025                     for (int k = 0; k < thisLookup.size(); k++) {  results[k].resize(thisLookup.size(), 0.0); }
1026
1027                     for (int j = 0; j < calcDists[i].size(); j++) {
1028                         int row = calcDists[i][j].seq1;
1029                         int column = calcDists[i][j].seq2;
1030                         double dist = calcDists[i][j].dist;
1031
1032                         results[row][column] = dist;
1033                         results[column][row] = dist;
1034                     }
1035                 }
1036             }
1037             for (int i = 0; i < calcDists.size(); i++) {  calcDists[i].clear(); }
1038                 }
1039
1040         if (iters != 0) {
1041             //we need to find the average distance and standard deviation for each groups distance
1042             vector< vector<seqDist>  > calcAverages = m->getAverages(calcDistsTotals, "average");
1043
1044             //print results
1045             for (int i = 0; i < calcDists.size(); i++) {
1046                 results.resize(thisLookup.size());
1047                 for (int k = 0; k < thisLookup.size(); k++) {  results[k].resize(thisLookup.size(), 0.0); }
1048
1049                 for (int j = 0; j < calcAverages[i].size(); j++) {
1050                     int row = calcAverages[i][j].seq1;
1051                     int column = calcAverages[i][j].seq2;
1052                     float dist = calcAverages[i][j].dist;
1053
1054                     results[row][column] = dist;
1055                     results[column][row] = dist;
1056                 }
1057             }
1058         }
1059
1060
1061         return results;
1062     }
1063     catch(exception& e) {
1064         m->errorOut(e, "GetMetaCommunityCommand", "generateDistanceMatrix");
1065         exit(1);
1066     }
1067 }
1068 /**************************************************************************************************/
1069 int GetMetaCommunityCommand::driver(vector<SharedRAbundVector*> thisLookup, vector< vector<seqDist> >& calcDists, Calculator* matrixCalculator) {
1070         try {
1071                 vector<SharedRAbundVector*> subset;
1072
1073                 for (int k = 0; k < thisLookup.size(); k++) { // pass cdd each set of groups to compare
1074
1075                         for (int l = 0; l < k; l++) {
1076
1077                                 if (k != l) { //we dont need to similiarity of a groups to itself
1078                                         subset.clear(); //clear out old pair of sharedrabunds
1079                                         //add new pair of sharedrabunds
1080                                         subset.push_back(thisLookup[k]); subset.push_back(thisLookup[l]);
1081
1082
1083
1084                     //if this calc needs all groups to calculate the pair load all groups
1085                     if (matrixCalculator->getNeedsAll()) {
1086                         //load subset with rest of lookup for those calcs that need everyone to calc for a pair
1087                         for (int w = 0; w < thisLookup.size(); w++) {
1088                             if ((w != k) && (w != l)) { subset.push_back(thisLookup[w]); }
1089                         }
1090                     }
1091
1092                     vector<double> tempdata = matrixCalculator->getValues(subset); //saves the calculator outputs
1093
1094                     if (m->control_pressed) { return 1; }
1095
1096                     seqDist temp(l, k, tempdata[0]);
1097                     //cout << l << '\t' << k << '\t' <<  tempdata[0] << endl;
1098                     calcDists[0].push_back(temp);
1099                 }
1100
1101                         }
1102                 }
1103
1104                 return 0;
1105         }
1106         catch(exception& e) {
1107                 m->errorOut(e, "MatrixOutputCommand", "driver");
1108                 exit(1);
1109         }
1110 }
1111 //**********************************************************************************************************************
1112