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             string relabund = relabunds[i];
 615             string matrixName = matrix[i];
 616             outputNames.push_back(matrixName); outputTypes["matrix"].push_back(matrixName);
 617
 618             finder->printZMatrix(matrixName, thisGroups);
 619
 620             double chi; vector<double> silhouettes;
 621             if (method == "dmm") {
 622                 double laplace = finder->getLaplace();
 623                 if(laplace < minLaplace){
 624                     minPartition = numPartitions;
 625                     minLaplace = laplace;
 626                 }
 627             }else {
 628                 chi = finder->calcCHIndex(dists);
 629                 silhouettes = finder->calcSilhouettes(dists);
 630                 if (chi > minLaplace) { //save partition with maximum ch index score
 631                     minPartition = numPartitions;
 632                     minLaplace = chi;
 633                     minSilhouettes = silhouettes;
 634                 }
 635             }
 636
 637             if (method == "dmm") {
 638                 finder->printFitData(cout, minLaplace);
 639                 finder->printFitData(fitData);
 640                 finder->printRelAbund(relabund, m->currentSharedBinLabels);
 641                 outputNames.push_back(relabund); outputTypes["relabund"].push_back(relabund);
 642             }else if ((method == "pam") || (method == "kmeans")) { //print silouettes and ch values
 643                 finder->printSilData(cout, chi, silhouettes);
 644                 finder->printSilData(silData, chi, silhouettes);
 645                 if (method == "kmeans") {
 646                     finder->printRelAbund(relabund, m->currentSharedBinLabels);
 647                     outputNames.push_back(relabund); outputTypes["relabund"].push_back(relabund);
 648                 }
 649             }
 650             delete finder;
 651
 652             if(optimizegap != -1 && (numPartitions - minPartition) >= optimizegap && numPartitions >= minpartitions){
 653                 string tempDoneFile = m->getRootName(m->getSimpleName(sharedfile)) + toString(processID) + ".done.temp";
 654                 ofstream outDone;
 655                 m->openOutputFile(tempDoneFile, outDone);
 656                 outDone << minPartition << endl;
 657                 outDone.close();
 658                 break;
 659             }
 660         }
 661         if (method == "dmm") { fitData.close(); }
 662
 663         if (m->control_pressed) { return 0; }
 664
 665         return minPartition;
 666     }
 667         catch(exception& e) {
 668                 m->errorOut(e, "GetMetaCommunityCommand", "processDriver");
 669                 exit(1);
 670         }
 671 }
 672 /**************************************************************************************************/
 673
 674 vector<double> GetMetaCommunityCommand::generateDesignFile(int numPartitions, map<string,string> variables){
 675     try {
 676         vector<double> piValues(numPartitions, 0);
 677
 678         ifstream postFile;
 679         variables["[tag]"] = toString(numPartitions);
 680         string input = getOutputFileName("matrix", variables);
 681         m->openInputFile(input, postFile);//((fileRoot + toString(numPartitions) + "mix.posterior").c_str()); //matrix file
 682         variables.erase("[tag]");
 683                 string outputFileName = getOutputFileName("design", variables);
 684         ofstream designFile;
 685         m->openOutputFile(outputFileName, designFile);
 686         outputNames.push_back(outputFileName); outputTypes["design"].push_back(outputFileName);
 687
 688
 689         vector<string> titles(numPartitions);
 690
 691         for(int i=0;i<numPartitions;i++){   postFile >> titles[i];  }
 692
 693         double posterior;
 694         string sampleName;
 695         int numSamples = 0;
 696
 697         while(postFile){
 698
 699             if (m->control_pressed) { break; }
 700
 701             double maxPosterior = 0.0000;
 702             int maxPartition = -1;
 703
 704             postFile >> sampleName;
 705
 706             for(int i=0;i<numPartitions;i++){
 707
 708                 postFile >> posterior;
 709                 if(posterior > maxPosterior){
 710                     maxPosterior = posterior;
 711                     maxPartition = i;
 712                 }
 713                 piValues[i] += posterior;
 714
 715             }
 716
 717             designFile << sampleName << '\t' << titles[maxPartition] << endl;
 718
 719             numSamples++;
 720             m->gobble(postFile);
 721         }
 722         for(int i=0;i<numPartitions;i++){
 723             piValues[i] /= (double)numSamples;
 724         }
 725
 726
 727         postFile.close();
 728         designFile.close();
 729
 730         return piValues;
 731     }
 732         catch(exception& e) {
 733                 m->errorOut(e, "GetMetaCommunityCommand", "generateDesignFile");
 734                 exit(1);
 735         }
 736 }
 737
 738 /**************************************************************************************************/
 739
 740 inline bool summaryFunction(summaryData i, summaryData j){ return i.difference > j.difference;   }
 741
 742 /**************************************************************************************************/
 743 int GetMetaCommunityCommand::generateSummaryFile(int numPartitions, map<string,string> v, vector<double> piValues){
 744     try {
 745         vector<summaryData> summary;
 746
 747         vector<double> pMean(numPartitions, 0);
 748         vector<double> pLCI(numPartitions, 0);
 749         vector<double> pUCI(numPartitions, 0);
 750
 751         string name, header;
 752         double mean, lci, uci;
 753
 754         ifstream referenceFile;
 755         map<string, string> variables;
 756         variables["[filename]"] = v["[filename]"];
 757         variables["[distance]"] = v["[distance]"];
 758         variables["[tag]"] = "1";
 759         string reference = getOutputFileName("relabund", variables);
 760         m->openInputFile(reference, referenceFile); //((fileRoot + label + ".1mix.relabund").c_str());
 761         variables["[tag]"] = toString(numPartitions);
 762         string partFile = getOutputFileName("relabund", variables);
 763         ifstream partitionFile;
 764         m->openInputFile(partFile, partitionFile); //((fileRoot + toString(numPartitions) + "mix.relabund").c_str());
 765
 766         header = m->getline(referenceFile);
 767         header = m->getline(partitionFile);
 768         stringstream head(header);
 769         string dummy, label;
 770         head >> dummy;
 771         vector<string> thetaValues(numPartitions, "");
 772         for(int i=0;i<numPartitions;i++){
 773             head >> label >> dummy >> dummy;
 774             thetaValues[i] = label.substr(label.find_last_of('_')+1);
 775         }
 776
 777
 778         vector<double> partitionDiff(numPartitions, 0.0000);
 779
 780         while(referenceFile){
 781
 782             if (m->control_pressed) { break; }
 783
 784             referenceFile >> name >> mean >> lci >> uci;
 785
 786             summaryData tempData;
 787             tempData.name = name;
 788             tempData.refMean = mean;
 789
 790             double difference = 0.0000;
 791
 792             partitionFile >> name;
 793             for(int j=0;j<numPartitions;j++){
 794                 partitionFile >> pMean[j] >> pLCI[j] >> pUCI[j];
 795                 difference += abs(mean - pMean[j]);
 796                 partitionDiff[j] += abs(mean - pMean[j]);;
 797             }
 798
 799             tempData.partMean = pMean;
 800             tempData.partLCI = pLCI;
 801             tempData.partUCI = pUCI;
 802             tempData.difference = difference;
 803             summary.push_back(tempData);
 804
 805             m->gobble(referenceFile);
 806             m->gobble(partitionFile);
 807         }
 808         referenceFile.close();
 809         partitionFile.close();
 810
 811         if (m->control_pressed) { return 0; }
 812
 813         int numOTUs = (int)summary.size();
 814
 815         sort(summary.begin(), summary.end(), summaryFunction);
 816
 817         variables.erase("[tag]");
 818                 string outputFileName = getOutputFileName("parameters", variables);
 819         outputNames.push_back(outputFileName); outputTypes["parameters"].push_back(outputFileName);
 820
 821         ofstream parameterFile;
 822         m->openOutputFile(outputFileName, parameterFile); //((fileRoot + "mix.parameters").c_str());
 823         parameterFile.setf(ios::fixed, ios::floatfield);
 824         parameterFile.setf(ios::showpoint);
 825
 826         double totalDifference =  0.0000;
 827         parameterFile << "Part\tDif2Ref_i\ttheta_i\tpi_i\n";
 828         for(int i=0;i<numPartitions;i++){
 829             if (m->control_pressed) { break; }
 830             parameterFile << i+1 << '\t' << setprecision(2) << partitionDiff[i] << '\t' << thetaValues[i] << '\t' << piValues[i] << endl;
 831             totalDifference += partitionDiff[i];
 832         }
 833         parameterFile.close();
 834
 835         if (m->control_pressed) { return 0; }
 836
 837         string summaryFileName = getOutputFileName("summary", variables);
 838         outputNames.push_back(summaryFileName); outputTypes["summary"].push_back(summaryFileName);
 839
 840         ofstream summaryFile;
 841         m->openOutputFile(summaryFileName, summaryFile); //((fileRoot + "mix.summary").c_str());
 842         summaryFile.setf(ios::fixed, ios::floatfield);
 843         summaryFile.setf(ios::showpoint);
 844
 845
 846         summaryFile << "OTU\tP0.mean";
 847         for(int i=0;i<numPartitions;i++){
 848             summaryFile << "\tP" << i+1 << ".mean\tP" << i+1 << ".lci\tP" << i+1 << ".uci";
 849         }
 850         summaryFile << "\tDifference\tCumFraction" << endl;
 851
 852         double cumDiff = 0.0000;
 853
 854         for(int i=0;i<numOTUs;i++){
 855             if (m->control_pressed) { break; }
 856             summaryFile << summary[i].name << setprecision(2) << '\t' << summary[i].refMean;
 857             for(int j=0;j<numPartitions;j++){
 858                 summaryFile  << '\t' << summary[i].partMean[j] << '\t' << summary[i].partLCI[j] << '\t' << summary[i].partUCI[j];
 859             }
 860
 861             cumDiff += summary[i].difference/totalDifference;
 862             summaryFile << '\t' << summary[i].difference << '\t' << cumDiff << endl;
 863         }
 864         summaryFile.close();
 865
 866         return 0;
 867
 868     }
 869         catch(exception& e) {
 870                 m->errorOut(e, "GetMetaCommunityCommand", "generateSummaryFile");
 871                 exit(1);
 872         }
 873
 874 }
 875 //**********************************************************************************************************************
 876 vector<vector<double> > GetMetaCommunityCommand::generateDistanceMatrix(vector<SharedRAbundVector*>& thisLookup){
 877     try {
 878         vector<vector<double> > results;
 879
 880         Calculator* matrixCalculator;
 881         ValidCalculators validCalculator;
 882         int i = 0;
 883
 884         if (validCalculator.isValidCalculator("matrix", Estimators[i]) == true) {
 885             if (Estimators[i] == "sharedsobs") {
 886                 matrixCalculator = new SharedSobsCS();
 887             }else if (Estimators[i] == "sharedchao") {
 888                 matrixCalculator = new SharedChao1();
 889             }else if (Estimators[i] == "sharedace") {
 890                 matrixCalculator = new SharedAce();
 891             }else if (Estimators[i] == "jabund") {
 892                 matrixCalculator = new JAbund();
 893             }else if (Estimators[i] == "sorabund") {
 894                 matrixCalculator = new SorAbund();
 895             }else if (Estimators[i] == "jclass") {
 896                 matrixCalculator = new Jclass();
 897             }else if (Estimators[i] == "sorclass") {
 898                 matrixCalculator = new SorClass();
 899             }else if (Estimators[i] == "jest") {
 900                 matrixCalculator = new Jest();
 901             }else if (Estimators[i] == "sorest") {
 902                 matrixCalculator = new SorEst();
 903             }else if (Estimators[i] == "thetayc") {
 904                 matrixCalculator = new ThetaYC();
 905             }else if (Estimators[i] == "thetan") {
 906                 matrixCalculator = new ThetaN();
 907             }else if (Estimators[i] == "kstest") {
 908                 matrixCalculator = new KSTest();
 909             }else if (Estimators[i] == "sharednseqs") {
 910                 matrixCalculator = new SharedNSeqs();
 911             }else if (Estimators[i] == "ochiai") {
 912                 matrixCalculator = new Ochiai();
 913             }else if (Estimators[i] == "anderberg") {
 914                 matrixCalculator = new Anderberg();
 915             }else if (Estimators[i] == "kulczynski") {
 916                 matrixCalculator = new Kulczynski();
 917             }else if (Estimators[i] == "kulczynskicody") {
 918                 matrixCalculator = new KulczynskiCody();
 919             }else if (Estimators[i] == "lennon") {
 920                 matrixCalculator = new Lennon();
 921             }else if (Estimators[i] == "morisitahorn") {
 922                 matrixCalculator = new MorHorn();
 923             }else if (Estimators[i] == "braycurtis") {
 924                 matrixCalculator = new BrayCurtis();
 925             }else if (Estimators[i] == "whittaker") {
 926                 matrixCalculator = new Whittaker();
 927             }else if (Estimators[i] == "odum") {
 928                 matrixCalculator = new Odum();
 929             }else if (Estimators[i] == "canberra") {
 930                 matrixCalculator = new Canberra();
 931             }else if (Estimators[i] == "structeuclidean") {
 932                 matrixCalculator = new StructEuclidean();
 933             }else if (Estimators[i] == "structchord") {
 934                 matrixCalculator = new StructChord();
 935             }else if (Estimators[i] == "hellinger") {
 936                 matrixCalculator = new Hellinger();
 937             }else if (Estimators[i] == "manhattan") {
 938                 matrixCalculator = new Manhattan();
 939             }else if (Estimators[i] == "structpearson") {
 940                 matrixCalculator = new StructPearson();
 941             }else if (Estimators[i] == "soergel") {
 942                 matrixCalculator = new Soergel();
 943             }else if (Estimators[i] == "spearman") {
 944                 matrixCalculator = new Spearman();
 945             }else if (Estimators[i] == "structkulczynski") {
 946                 matrixCalculator = new StructKulczynski();
 947             }else if (Estimators[i] == "speciesprofile") {
 948                 matrixCalculator = new SpeciesProfile();
 949             }else if (Estimators[i] == "hamming") {
 950                 matrixCalculator = new Hamming();
 951             }else if (Estimators[i] == "structchi2") {
 952                 matrixCalculator = new StructChi2();
 953             }else if (Estimators[i] == "gower") {
 954                 matrixCalculator = new Gower();
 955             }else if (Estimators[i] == "memchi2") {
 956                 matrixCalculator = new MemChi2();
 957             }else if (Estimators[i] == "memchord") {
 958                 matrixCalculator = new MemChord();
 959             }else if (Estimators[i] == "memeuclidean") {
 960                 matrixCalculator = new MemEuclidean();
 961             }else if (Estimators[i] == "mempearson") {
 962                 matrixCalculator = new MemPearson();
 963             }else if (Estimators[i] == "jsd") {
 964                 matrixCalculator = new JSD();
 965             }else if (Estimators[i] == "rjsd") {
 966                 matrixCalculator = new RJSD();
 967             }else {
 968                 m->mothurOut("[ERROR]: " + Estimators[i] + " is not a valid calculator, please correct.\n"); m->control_pressed = true; return results;
 969             }
 970         }
 971
 972         //calc distances
 973         vector< vector< vector<seqDist> > > calcDistsTotals;  //each iter, then each groupCombos dists. this will be used to make .dist files
 974         vector< vector<seqDist> > calcDists; calcDists.resize(1);
 975
 976         for (int thisIter = 0; thisIter < iters+1; thisIter++) {
 977
 978             vector<SharedRAbundVector*> thisItersLookup = thisLookup;
 979
 980             if (subsample && (thisIter != 0)) {
 981                 SubSample sample;
 982                 vector<string> tempLabels; //dont need since we arent printing the sampled sharedRabunds
 983
 984                 //make copy of lookup so we don't get access violations
 985                 vector<SharedRAbundVector*> newLookup;
 986                 for (int k = 0; k < thisItersLookup.size(); k++) {
 987                     SharedRAbundVector* temp = new SharedRAbundVector();
 988                     temp->setLabel(thisItersLookup[k]->getLabel());
 989                     temp->setGroup(thisItersLookup[k]->getGroup());
 990                     newLookup.push_back(temp);
 991                 }
 992
 993                 //for each bin
 994                 for (int k = 0; k < thisItersLookup[0]->getNumBins(); k++) {
 995                     if (m->control_pressed) { for (int j = 0; j < newLookup.size(); j++) {  delete newLookup[j];  } return results; }
 996                     for (int j = 0; j < thisItersLookup.size(); j++) { newLookup[j]->push_back(thisItersLookup[j]->getAbundance(k), thisItersLookup[j]->getGroup()); }
 997                 }
 998
 999                 tempLabels = sample.getSample(newLookup, subsampleSize);
1000                 thisItersLookup = newLookup;
1001             }
1002
1003
1004             driver(thisItersLookup, calcDists, matrixCalculator);
1005
1006             if (subsample && (thisIter != 0)) {
1007                 if((thisIter) % 100 == 0){      m->mothurOutJustToScreen(toString(thisIter)+"\n");              }
1008                 calcDistsTotals.push_back(calcDists);
1009                 for (int i = 0; i < calcDists.size(); i++) {
1010                     for (int j = 0; j < calcDists[i].size(); j++) {
1011                         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");  }
1012                     }
1013                 }
1014                 //clean up memory
1015                 for (int i = 0; i < thisItersLookup.size(); i++) { delete thisItersLookup[i]; }
1016                 thisItersLookup.clear();
1017             }else { //print results for whole dataset
1018                 for (int i = 0; i < calcDists.size(); i++) {
1019                     if (m->control_pressed) { break; }
1020
1021                     //initialize matrix
1022                     results.resize(thisLookup.size());
1023                     for (int k = 0; k < thisLookup.size(); k++) {  results[k].resize(thisLookup.size(), 0.0); }
1024
1025                     for (int j = 0; j < calcDists[i].size(); j++) {
1026                         int row = calcDists[i][j].seq1;
1027                         int column = calcDists[i][j].seq2;
1028                         double dist = calcDists[i][j].dist;
1029
1030                         results[row][column] = dist;
1031                         results[column][row] = dist;
1032                     }
1033                 }
1034             }
1035             for (int i = 0; i < calcDists.size(); i++) {  calcDists[i].clear(); }
1036                 }
1037
1038         if (iters != 0) {
1039             //we need to find the average distance and standard deviation for each groups distance
1040             vector< vector<seqDist>  > calcAverages = m->getAverages(calcDistsTotals, "average");
1041
1042             //print results
1043             for (int i = 0; i < calcDists.size(); i++) {
1044                 results.resize(thisLookup.size());
1045                 for (int k = 0; k < thisLookup.size(); k++) {  results[k].resize(thisLookup.size(), 0.0); }
1046
1047                 for (int j = 0; j < calcAverages[i].size(); j++) {
1048                     int row = calcAverages[i][j].seq1;
1049                     int column = calcAverages[i][j].seq2;
1050                     float dist = calcAverages[i][j].dist;
1051
1052                     results[row][column] = dist;
1053                     results[column][row] = dist;
1054                 }
1055             }
1056         }
1057
1058
1059         return results;
1060     }
1061     catch(exception& e) {
1062         m->errorOut(e, "GetMetaCommunityCommand", "generateDistanceMatrix");
1063         exit(1);
1064     }
1065 }
1066 /**************************************************************************************************/
1067 int GetMetaCommunityCommand::driver(vector<SharedRAbundVector*> thisLookup, vector< vector<seqDist> >& calcDists, Calculator* matrixCalculator) {
1068         try {
1069                 vector<SharedRAbundVector*> subset;
1070
1071                 for (int k = 0; k < thisLookup.size(); k++) { // pass cdd each set of groups to compare
1072
1073                         for (int l = 0; l < k; l++) {
1074
1075                                 if (k != l) { //we dont need to similiarity of a groups to itself
1076                                         subset.clear(); //clear out old pair of sharedrabunds
1077                                         //add new pair of sharedrabunds
1078                                         subset.push_back(thisLookup[k]); subset.push_back(thisLookup[l]);
1079
1080
1081
1082                     //if this calc needs all groups to calculate the pair load all groups
1083                     if (matrixCalculator->getNeedsAll()) {
1084                         //load subset with rest of lookup for those calcs that need everyone to calc for a pair
1085                         for (int w = 0; w < thisLookup.size(); w++) {
1086                             if ((w != k) && (w != l)) { subset.push_back(thisLookup[w]); }
1087                         }
1088                     }
1089
1090                     vector<double> tempdata = matrixCalculator->getValues(subset); //saves the calculator outputs
1091
1092                     if (m->control_pressed) { return 1; }
1093
1094                     seqDist temp(l, k, tempdata[0]);
1095                     //cout << l << '\t' << k << '\t' <<  tempdata[0] << endl;
1096                     calcDists[0].push_back(temp);
1097                 }
1098
1099                         }
1100                 }
1101
1102                 return 0;
1103         }
1104         catch(exception& e) {
1105                 m->errorOut(e, "MatrixOutputCommand", "driver");
1106                 exit(1);
1107         }
1108 }
1109 //**********************************************************************************************************************
1110