amovacommand.cpp

   1 /*
   2  *  amovacommand.cpp
   3  *  mothur
   4  *
   5  *  Created by westcott on 2/7/11.
   6  *  Copyright 2011 Schloss Lab. All rights reserved.
   7  *
   8  */
   9
  10 #include "amovacommand.h"
  11 #include "sharedutilities.h"
  12 #include "sharedsobscollectsummary.h"
  13 #include "sharedchao1.h"
  14 #include "sharedace.h"
  15 #include "sharednseqs.h"
  16 #include "sharedjabund.h"
  17 #include "sharedsorabund.h"
  18 #include "sharedjclass.h"
  19 #include "sharedsorclass.h"
  20 #include "sharedjest.h"
  21 #include "sharedsorest.h"
  22 #include "sharedthetayc.h"
  23 #include "sharedthetan.h"
  24 #include "sharedkstest.h"
  25 #include "whittaker.h"
  26 #include "sharedochiai.h"
  27 #include "sharedanderbergs.h"
  28 #include "sharedkulczynski.h"
  29 #include "sharedkulczynskicody.h"
  30 #include "sharedlennon.h"
  31 #include "sharedmorisitahorn.h"
  32 #include "sharedbraycurtis.h"
  33 #include "sharedjackknife.h"
  34 #include "whittaker.h"
  35 #include "odum.h"
  36 #include "canberra.h"
  37 #include "structeuclidean.h"
  38 #include "structchord.h"
  39 #include "hellinger.h"
  40 #include "manhattan.h"
  41 #include "structpearson.h"
  42 #include "soergel.h"
  43 #include "spearman.h"
  44 #include "structkulczynski.h"
  45 #include "structchi2.h"
  46 #include "speciesprofile.h"
  47 #include "hamming.h"
  48 #include "gower.h"
  49 #include "memchi2.h"
  50 #include "memchord.h"
  51 #include "memeuclidean.h"
  52 #include "mempearson.h"
  53
  54 //**********************************************************************************************************************
  55 vector<string> AmovaCommand::getValidParameters(){
  56         try {
  57                 string Array[] =  {"groups","label","outputdir","iters","phylip","design","sets","processors","inputdir"};
  58                 vector<string> myArray (Array, Array+(sizeof(Array)/sizeof(string)));
  59                 return myArray;
  60         }
  61         catch(exception& e) {
  62                 m->errorOut(e, "AmovaCommand", "getValidParameters");
  63                 exit(1);
  64         }
  65 }
  66 //**********************************************************************************************************************
  67 AmovaCommand::AmovaCommand(){
  68         try {
  69                 abort = true; calledHelp = true;
  70                 vector<string> tempOutNames;
  71                 outputTypes["amova"] = tempOutNames;
  72         }
  73         catch(exception& e) {
  74                 m->errorOut(e, "AmovaCommand", "AmovaCommand");
  75                 exit(1);
  76         }
  77 }
  78 //**********************************************************************************************************************
  79 vector<string> AmovaCommand::getRequiredParameters(){
  80         try {
  81                 string Array[] =  {"design"};
  82                 vector<string> myArray (Array, Array+(sizeof(Array)/sizeof(string)));
  83                 return myArray;
  84         }
  85         catch(exception& e) {
  86                 m->errorOut(e, "AmovaCommand", "getRequiredParameters");
  87                 exit(1);
  88         }
  89 }
  90 //**********************************************************************************************************************
  91 vector<string> AmovaCommand::getRequiredFiles(){
  92         try {
  93                 string Array[] =  {};
  94                 vector<string> myArray (Array, Array+(sizeof(Array)/sizeof(string)));
  95                 return myArray;
  96         }
  97         catch(exception& e) {
  98                 m->errorOut(e, "AmovaCommand", "getRequiredFiles");
  99                 exit(1);
 100         }
 101 }
 102 //**********************************************************************************************************************
 103
 104 AmovaCommand::AmovaCommand(string option) {
 105         try {
 106                 globaldata = GlobalData::getInstance();
 107                 abort = false; calledHelp = false;
 108                 allLines = 1;
 109                 labels.clear();
 110
 111                 //allow user to run help
 112                 if(option == "help") { help(); abort = true; calledHelp = true; }
 113
 114                 else {
 115                         //valid paramters for this command
 116                         string AlignArray[] =  {"groups","label","outputdir","iters","phylip","design","sets","processors","inputdir"};
 117                         vector<string> myArray (AlignArray, AlignArray+(sizeof(AlignArray)/sizeof(string)));
 118
 119                         OptionParser parser(option);
 120                         map<string,string> parameters = parser.getParameters();
 121
 122                         ValidParameters validParameter;
 123
 124                         //check to make sure all parameters are valid for command
 125                         map<string,string>::iterator it;
 126                         for (it = parameters.begin(); it != parameters.end(); it++) {
 127                                 if (validParameter.isValidParameter(it->first, myArray, it->second) != true) {  abort = true;  }
 128                         }
 129
 130                         //initialize outputTypes
 131                         vector<string> tempOutNames;
 132                         outputTypes["amova"] = tempOutNames;
 133
 134                         //if the user changes the output directory command factory will send this info to us in the output parameter
 135                         outputDir = validParameter.validFile(parameters, "outputdir", false);           if (outputDir == "not found"){  outputDir = ""; }
 136
 137                         //if the user changes the input directory command factory will send this info to us in the output parameter
 138                         string inputDir = validParameter.validFile(parameters, "inputdir", false);
 139                         if (inputDir == "not found"){   inputDir = "";          }
 140                         else {
 141                                 string path;
 142                                 it = parameters.find("design");
 143                                 //user has given a template file
 144                                 if(it != parameters.end()){
 145                                         path = m->hasPath(it->second);
 146                                         //if the user has not given a path then, add inputdir. else leave path alone.
 147                                         if (path == "") {       parameters["design"] = inputDir + it->second;           }
 148                                 }
 149
 150                                 it = parameters.find("phylip");
 151                                 //user has given a template file
 152                                 if(it != parameters.end()){
 153                                         path = m->hasPath(it->second);
 154                                         //if the user has not given a path then, add inputdir. else leave path alone.
 155                                         if (path == "") {       parameters["phylip"] = inputDir + it->second;           }
 156                                 }
 157                         }
 158
 159                         phylipfile = validParameter.validFile(parameters, "phylip", true);
 160                         if (phylipfile == "not open") { phylipfile = ""; abort = true; }
 161                         else if (phylipfile == "not found") { phylipfile = ""; }
 162                         else {  globaldata->newRead(); format = "phylip";  globaldata->setPhylipFile(phylipfile);       }
 163
 164                         //check for required parameters
 165                         designfile = validParameter.validFile(parameters, "design", true);
 166                         if (designfile == "not open") { abort = true; }
 167                         else if (designfile == "not found") {  designfile = "";  m->mothurOut("You must provide an design file."); m->mothurOutEndLine(); abort = true; }
 168
 169                         //make sure the user has already run the read.otu command
 170                         if ((globaldata->getSharedFile() == "")) {
 171                                 if ((phylipfile == "")) {
 172                                         m->mothurOut("You must read a list and a group, a shared file or provide a distance matrix before you can use the amova command."); m->mothurOutEndLine(); abort = true;
 173                                 }
 174                         }else { sharedfile = globaldata->getSharedFile(); }
 175
 176                         //use distance matrix if one is provided
 177                         if ((sharedfile != "") && (phylipfile != "")) { sharedfile = ""; }
 178
 179                         //check for optional parameter and set defaults
 180                         // ...at some point should added some additional type checking...
 181                         label = validParameter.validFile(parameters, "label", false);
 182                         if (label == "not found") { label = ""; }
 183                         else {
 184                                 if(label != "all") {  m->splitAtDash(label, labels);  allLines = 0;  }
 185                                 else { allLines = 1;  }
 186                         }
 187
 188                         //if the user has not specified any labels use the ones from read.otu
 189                         if (label == "") {
 190                                 allLines = globaldata->allLines;
 191                                 labels = globaldata->labels;
 192                         }
 193
 194                         groups = validParameter.validFile(parameters, "groups", false);
 195                         if (groups == "not found") { groups = "";  }
 196                         else {
 197                                 m->splitAtDash(groups, Groups);
 198                                 globaldata->Groups = Groups;
 199                         }
 200
 201                         sets = validParameter.validFile(parameters, "sets", false);
 202                         if (sets == "not found") { sets = ""; pickedGroups = false; }
 203                         else {
 204                                 pickedGroups = true;
 205                                 m->splitAtDash(sets, Sets);
 206                         }
 207
 208
 209                         string temp = validParameter.validFile(parameters, "iters", false);                     if (temp == "not found") { temp = "1000"; }
 210                         convert(temp, iters);
 211
 212                         temp = validParameter.validFile(parameters, "processors", false);                       if (temp == "not found"){       temp = "1";     }
 213                         convert(temp, processors);
 214
 215                         vector<string> Estimators;
 216                         calc = validParameter.validFile(parameters, "calc", false);
 217                         if (calc == "not found") { calc = "morisitahorn";  }
 218                         m->splitAtDash(calc, Estimators);
 219
 220                         if (abort == false) {
 221
 222                                 ValidCalculators* validCalculator = new ValidCalculators();
 223
 224                                 for (int i=0; i<Estimators.size(); i++) {
 225                                         if (validCalculator->isValidCalculator("treegroup", Estimators[i]) == true) {
 226                                                 if (Estimators[i] == "sharedsobs") {
 227                                                         calculators.push_back(new SharedSobsCS());
 228                                                 }else if (Estimators[i] == "sharedchao") {
 229                                                         calculators.push_back(new SharedChao1());
 230                                                 }else if (Estimators[i] == "sharedace") {
 231                                                         calculators.push_back(new SharedAce());
 232                                                 }else if (Estimators[i] == "jabund") {
 233                                                         calculators.push_back(new JAbund());
 234                                                 }else if (Estimators[i] == "sorabund") {
 235                                                         calculators.push_back(new SorAbund());
 236                                                 }else if (Estimators[i] == "jclass") {
 237                                                         calculators.push_back(new Jclass());
 238                                                 }else if (Estimators[i] == "sorclass") {
 239                                                         calculators.push_back(new SorClass());
 240                                                 }else if (Estimators[i] == "jest") {
 241                                                         calculators.push_back(new Jest());
 242                                                 }else if (Estimators[i] == "sorest") {
 243                                                         calculators.push_back(new SorEst());
 244                                                 }else if (Estimators[i] == "thetayc") {
 245                                                         calculators.push_back(new ThetaYC());
 246                                                 }else if (Estimators[i] == "thetan") {
 247                                                         calculators.push_back(new ThetaN());
 248                                                 }else if (Estimators[i] == "kstest") {
 249                                                         calculators.push_back(new KSTest());
 250                                                 }else if (Estimators[i] == "sharednseqs") {
 251                                                         calculators.push_back(new SharedNSeqs());
 252                                                 }else if (Estimators[i] == "ochiai") {
 253                                                         calculators.push_back(new Ochiai());
 254                                                 }else if (Estimators[i] == "anderberg") {
 255                                                         calculators.push_back(new Anderberg());
 256                                                 }else if (Estimators[i] == "kulczynski") {
 257                                                         calculators.push_back(new Kulczynski());
 258                                                 }else if (Estimators[i] == "kulczynskicody") {
 259                                                         calculators.push_back(new KulczynskiCody());
 260                                                 }else if (Estimators[i] == "lennon") {
 261                                                         calculators.push_back(new Lennon());
 262                                                 }else if (Estimators[i] == "morisitahorn") {
 263                                                         calculators.push_back(new MorHorn());
 264                                                 }else if (Estimators[i] == "braycurtis") {
 265                                                         calculators.push_back(new BrayCurtis());
 266                                                 }else if (Estimators[i] == "whittaker") {
 267                                                         calculators.push_back(new Whittaker());
 268                                                 }else if (Estimators[i] == "odum") {
 269                                                         calculators.push_back(new Odum());
 270                                                 }else if (Estimators[i] == "canberra") {
 271                                                         calculators.push_back(new Canberra());
 272                                                 }else if (Estimators[i] == "structeuclidean") {
 273                                                         calculators.push_back(new StructEuclidean());
 274                                                 }else if (Estimators[i] == "structchord") {
 275                                                         calculators.push_back(new StructChord());
 276                                                 }else if (Estimators[i] == "hellinger") {
 277                                                         calculators.push_back(new Hellinger());
 278                                                 }else if (Estimators[i] == "manhattan") {
 279                                                         calculators.push_back(new Manhattan());
 280                                                 }else if (Estimators[i] == "structpearson") {
 281                                                         calculators.push_back(new StructPearson());
 282                                                 }else if (Estimators[i] == "soergel") {
 283                                                         calculators.push_back(new Soergel());
 284                                                 }else if (Estimators[i] == "spearman") {
 285                                                         calculators.push_back(new Spearman());
 286                                                 }else if (Estimators[i] == "structkulczynski") {
 287                                                         calculators.push_back(new StructKulczynski());
 288                                                 }else if (Estimators[i] == "speciesprofile") {
 289                                                         calculators.push_back(new SpeciesProfile());
 290                                                 }else if (Estimators[i] == "hamming") {
 291                                                         calculators.push_back(new Hamming());
 292                                                 }else if (Estimators[i] == "structchi2") {
 293                                                         calculators.push_back(new StructChi2());
 294                                                 }else if (Estimators[i] == "gower") {
 295                                                         calculators.push_back(new Gower());
 296                                                 }else if (Estimators[i] == "memchi2") {
 297                                                         calculators.push_back(new MemChi2());
 298                                                 }else if (Estimators[i] == "memchord") {
 299                                                         calculators.push_back(new MemChord());
 300                                                 }else if (Estimators[i] == "memeuclidean") {
 301                                                         calculators.push_back(new MemEuclidean());
 302                                                 }else if (Estimators[i] == "mempearson") {
 303                                                         calculators.push_back(new MemPearson());
 304                                                 }
 305                                         }
 306                                 }
 307                         }
 308                 }
 309
 310         }
 311         catch(exception& e) {
 312                 m->errorOut(e, "AmovaCommand", "AmovaCommand");
 313                 exit(1);
 314         }
 315 }
 316
 317 //**********************************************************************************************************************
 318
 319 void AmovaCommand::help(){
 320         try {
 321                 m->mothurOut("The amova command can only be executed after a successful read.otu command of a list and group or shared file, or by providing a phylip formatted distance matrix.\n");
 322                 m->mothurOut("The amova command outputs a .amova file. \n");
 323                 m->mothurOut("The amova command parameters are phylip, iters, groups, label, design, sets and processors.  The design parameter is required.\n");
 324                 m->mothurOut("The design parameter allows you to assign your groups to sets when you are running amova. It is required. \n");
 325                 m->mothurOut("The design file looks like the group file.  It is a 2 column tab delimited file, where the first column is the group name and the second column is the set the group belongs to.\n");
 326                 m->mothurOut("The sets parameter allows you to specify which of the sets in your designfile you would like to analyze. The set names are separated by dashes. THe default is all sets in the designfile. To run the pairwise comparisons of all sets use sets=all.\n");
 327                 m->mothurOut("The iters parameter allows you to set number of randomization for the P value.  The default is 1000. \n");
 328                 m->mothurOut("The groups parameter allows you to specify which of the groups you would like included. The group names are separated by dashes. groups=all will find all pairwise comparisons. \n");
 329                 m->mothurOut("The label parameter allows you to select what distance levels you would like, and are also separated by dashes.\n");
 330                 m->mothurOut("The processors parameter allows you to specify how many processors you would like to use.  The default is 1. \n");
 331                 m->mothurOut("The amova command should be in the following format: amova(design=yourDesignFile).\n");
 332                 m->mothurOut("Example amova(design=temp.design, groups=A-B-C).\n");
 333                 m->mothurOut("Note: No spaces between parameter labels (i.e. groups), '=' and parameters (i.e.yourGroups).\n\n");
 334
 335         }
 336         catch(exception& e) {
 337                 m->errorOut(e, "AmovaCommand", "help");
 338                 exit(1);
 339         }
 340 }
 341
 342 //**********************************************************************************************************************
 343
 344 AmovaCommand::~AmovaCommand(){}
 345
 346 //**********************************************************************************************************************
 347
 348 int AmovaCommand::execute(){
 349         try {
 350
 351                 if (abort == true) { if (calledHelp) { return 0; }  return 2;   }
 352
 353                 //read design file
 354                 designMap = new GroupMap(designfile);
 355                 designMap->readDesignMap();
 356
 357                 //make sure sets are all in designMap
 358                 SharedUtil* util = new SharedUtil();
 359                 util->setGroups(Sets, designMap->namesOfGroups);
 360                 delete util;
 361
 362                 //if the user pickedGroups or set sets=all, then we want to figure out how to split the pairwise comparison between processors
 363                 int numGroups = Sets.size();
 364                 if (pickedGroups) {
 365                         for (int a=0; a<numGroups; a++) {
 366                                 for (int l = 0; l < a; l++) {
 367                                         vector<string> groups; groups.push_back(Sets[a]); groups.push_back(Sets[l]);
 368                                         namesOfGroupCombos.push_back(groups);
 369                                 }
 370                         }
 371                 }else { //one giant compare
 372                         vector<string> groups;
 373                         for (int a=0; a<numGroups; a++) { groups.push_back(Sets[a]); }
 374                         namesOfGroupCombos.push_back(groups);
 375                 }
 376
 377                 //only 1 combo
 378                 if (numGroups == 2) { processors = 1; }
 379                 else if (numGroups < 2) { m->mothurOut("Not enough sets, I need at least 2 valid sets. Unable to complete command."); m->mothurOutEndLine(); m->control_pressed = true; }
 380
 381                 #if defined (__APPLE__) || (__MACH__) || (linux) || (__linux)
 382                         if(processors != 1){
 383                                 int numPairs = namesOfGroupCombos.size();
 384                                 int numPairsPerProcessor = numPairs / processors;
 385
 386                                 for (int i = 0; i < processors; i++) {
 387                                         int startPos = i * numPairsPerProcessor;
 388                                         if(i == processors - 1){
 389                                                 numPairsPerProcessor = numPairs - i * numPairsPerProcessor;
 390                                         }
 391                                         lines.push_back(linePair(startPos, numPairsPerProcessor));
 392                                 }
 393                         }
 394                 #endif
 395
 396                 if (sharedfile != "") { //create distance matrix for each label
 397
 398                         if (outputDir == "") { outputDir = m->hasPath(sharedfile); }
 399
 400                         //for each calculator
 401                         for(int i = 0 ; i < calculators.size(); i++) {
 402
 403                                 //create a new filename
 404                                 ofstream out;
 405                                 string outputFile = outputDir + m->getRootName(m->getSimpleName(sharedfile)) + calculators[i]->getName() + ".amova";
 406                                 m->openOutputFile(outputFile, out);
 407                                 outputNames.push_back(outputFile); outputTypes["amova"].push_back(outputFile);
 408
 409                                 //print headers
 410                                 out << "label\tgroupsCompared\tMeanSquaredWithin\tMeanSquaredAmong\tFstatistic\tpValue" << endl;
 411                                 out.close();
 412                         }
 413
 414                         InputData input(sharedfile, "sharedfile");
 415                         vector<SharedRAbundVector*> lookup = input.getSharedRAbundVectors();
 416                         string lastLabel = lookup[0]->getLabel();
 417
 418                         //if the users enters label "0.06" and there is no "0.06" in their file use the next lowest label.
 419                         set<string> processedLabels;
 420                         set<string> userLabels = labels;
 421
 422                         //sanity check between shared file groups and design file groups
 423                         for (int i = 0; i < lookup.size(); i++) {
 424                                 string group = designMap->getGroup(lookup[i]->getGroup());
 425
 426                                 if (group == "not found") { m->control_pressed = true;  m->mothurOut("[ERROR]: " + lookup[i]->getGroup() + " is not in your design file, please correct."); m->mothurOutEndLine();  }
 427                         }
 428
 429                         //as long as you are not at the end of the file or done wih the lines you want
 430                         while((lookup[0] != NULL) && ((allLines == 1) || (userLabels.size() != 0))) {
 431
 432                                 if (m->control_pressed) {  for (int i = 0; i < lookup.size(); i++) {  delete lookup[i];  } globaldata->Groups.clear();  delete designMap;  for (int i = 0; i < outputNames.size(); i++) {       remove(outputNames[i].c_str()); } return 0; }
 433
 434                                 if(allLines == 1 || labels.count(lookup[0]->getLabel()) == 1){
 435
 436                                         m->mothurOut(lookup[0]->getLabel()); m->mothurOutEndLine();
 437                                         process(lookup);
 438
 439                                         processedLabels.insert(lookup[0]->getLabel());
 440                                         userLabels.erase(lookup[0]->getLabel());
 441                                 }
 442
 443                                 if ((m->anyLabelsToProcess(lookup[0]->getLabel(), userLabels, "") == true) && (processedLabels.count(lastLabel) != 1)) {
 444                                         string saveLabel = lookup[0]->getLabel();
 445
 446                                         for (int i = 0; i < lookup.size(); i++) {  delete lookup[i];  }
 447                                         lookup = input.getSharedRAbundVectors(lastLabel);
 448                                         m->mothurOut(lookup[0]->getLabel()); m->mothurOutEndLine();
 449
 450                                         process(lookup);
 451
 452                                         processedLabels.insert(lookup[0]->getLabel());
 453                                         userLabels.erase(lookup[0]->getLabel());
 454
 455                                         //restore real lastlabel to save below
 456                                         lookup[0]->setLabel(saveLabel);
 457                                 }
 458
 459                                 lastLabel = lookup[0]->getLabel();
 460                                 //prevent memory leak
 461                                 for (int i = 0; i < lookup.size(); i++) {  delete lookup[i]; lookup[i] = NULL; }
 462
 463                                 if (m->control_pressed) {  globaldata->Groups.clear();   delete designMap;  for (int i = 0; i < outputNames.size(); i++) {      remove(outputNames[i].c_str()); } return 0; }
 464
 465                                 //get next line to process
 466                                 lookup = input.getSharedRAbundVectors();
 467                         }
 468
 469                         if (m->control_pressed) { globaldata->Groups.clear();  delete designMap;  for (int i = 0; i < outputNames.size(); i++) {        remove(outputNames[i].c_str()); }  return 0; }
 470
 471                         //output error messages about any remaining user labels
 472                         set<string>::iterator it;
 473                         bool needToRun = false;
 474                         for (it = userLabels.begin(); it != userLabels.end(); it++) {
 475                                 m->mothurOut("Your file does not include the label " + *it);
 476                                 if (processedLabels.count(lastLabel) != 1) {
 477                                         m->mothurOut(". I will use " + lastLabel + "."); m->mothurOutEndLine();
 478                                         needToRun = true;
 479                                 }else {
 480                                         m->mothurOut(". Please refer to " + lastLabel + "."); m->mothurOutEndLine();
 481                                 }
 482                         }
 483
 484                         //run last label if you need to
 485                         if (needToRun == true)  {
 486                                 for (int i = 0; i < lookup.size(); i++) { if (lookup[i] != NULL) { delete lookup[i]; } }
 487                                 lookup = input.getSharedRAbundVectors(lastLabel);
 488
 489                                 m->mothurOut(lookup[0]->getLabel()); m->mothurOutEndLine();
 490
 491                                 process(lookup);
 492
 493                                 for (int i = 0; i < lookup.size(); i++) {  delete lookup[i];  }
 494                         }
 495
 496                         //reset groups parameter
 497                         globaldata->Groups.clear();
 498
 499
 500                 }else { //user provided distance matrix
 501
 502                         if (outputDir == "") { outputDir = m->hasPath(phylipfile); }
 503
 504                         //create a new filename
 505                         ofstream out;
 506                         string outputFile = outputDir + m->getRootName(m->getSimpleName(phylipfile))  + "amova";
 507                         m->openOutputFile(outputFile, out);
 508                         outputNames.push_back(outputFile); outputTypes["amova"].push_back(outputFile);
 509
 510                         //print headers
 511                         out << "groupsCompared\tMeanSquaredWithin\tMeanSquaredAmong\tFstatistic\tpValue" << endl;
 512                         out.close();
 513
 514                         ReadPhylipVector readMatrix(phylipfile);
 515                         vector< vector<double> > completeMatrix;
 516                         vector<string> names = readMatrix.read(completeMatrix);
 517
 518                         #if defined (__APPLE__) || (__MACH__) || (linux) || (__linux)
 519                                 if(processors == 1){
 520                                         driver(0, namesOfGroupCombos.size(), names, "", completeMatrix);
 521                                 }else{
 522                                         int process = 1;
 523                                         vector<int> processIDS;
 524
 525                                         //loop through and create all the processes you want
 526                                         while (process != processors) {
 527                                                 int pid = fork();
 528
 529                                                 if (pid > 0) {
 530                                                         processIDS.push_back(pid);
 531                                                         process++;
 532                                                 }else if (pid == 0){
 533                                                         driver(lines[process].start, lines[process].num, names, toString(getpid()), completeMatrix);
 534                                                         exit(0);
 535                                                 }else {
 536                                                         m->mothurOut("[ERROR]: unable to spawn the necessary processes."); m->mothurOutEndLine();
 537                                                         for (int i = 0; i < processIDS.size(); i++) { kill (processIDS[i], SIGINT); }
 538                                                         exit(0);
 539                                                 }
 540                                         }
 541
 542                                         //do my part
 543                                         driver(lines[0].start, lines[0].num, names, "", completeMatrix);
 544
 545                                         //force parent to wait until all the processes are done
 546                                         for (int i=0;i<(processors-1);i++) {
 547                                                 int temp = processIDS[i];
 548                                                 wait(&temp);
 549                                         }
 550
 551                                         //append files
 552                                         string outputFile = outputDir + m->getRootName(m->getSimpleName(phylipfile)) + "amova";
 553                                         for (int j = 0; j < processIDS.size(); j++) {
 554                                                 m->appendFiles((outputFile + toString(processIDS[j])), outputFile);
 555                                                 remove((outputFile + toString(processIDS[j])).c_str());
 556                                         }
 557
 558                                 }
 559                         #else
 560                                 driver(0, namesOfGroupCombos.size(), names, "", completeMatrix);
 561                         #endif
 562
 563                 }
 564
 565                 delete designMap;
 566
 567                 if (m->control_pressed) { for (int i = 0; i < outputNames.size(); i++) {        remove(outputNames[i].c_str()); } return 0; }
 568
 569                 m->mothurOutEndLine();
 570                 m->mothurOut("Output File Names: "); m->mothurOutEndLine();
 571                 for (int i = 0; i < outputNames.size(); i++) {  m->mothurOut(outputNames[i]); m->mothurOutEndLine();    }
 572                 m->mothurOutEndLine();
 573
 574                 return 0;
 575         }
 576         catch(exception& e) {
 577                 m->errorOut(e, "AmovaCommand", "execute");
 578                 exit(1);
 579         }
 580 }
 581 //**********************************************************************************************************************
 582
 583 int AmovaCommand::process(vector<SharedRAbundVector*> thisLookup) {
 584         try{
 585
 586 #if defined (__APPLE__) || (__MACH__) || (linux) || (__linux)
 587                 if(processors == 1){
 588                         driver(0, namesOfGroupCombos.size(), thisLookup, "");
 589                 }else{
 590                         int process = 1;
 591                         vector<int> processIDS;
 592
 593                         //loop through and create all the processes you want
 594                         while (process != processors) {
 595                                 int pid = fork();
 596
 597                                 if (pid > 0) {
 598                                         processIDS.push_back(pid);
 599                                         process++;
 600                                 }else if (pid == 0){
 601                                         driver(lines[process].start, lines[process].num, thisLookup, toString(getpid()));
 602                                         exit(0);
 603                                 }else {
 604                                         m->mothurOut("[ERROR]: unable to spawn the necessary processes."); m->mothurOutEndLine();
 605                                         for (int i = 0; i < processIDS.size(); i++) { kill (processIDS[i], SIGINT); }
 606                                         exit(0);
 607                                 }
 608                         }
 609
 610                         //do my part
 611                         driver(lines[0].start, lines[0].num, thisLookup, "");
 612
 613                         //force parent to wait until all the processes are done
 614                         for (int i=0;i<(processors-1);i++) {
 615                                 int temp = processIDS[i];
 616                                 wait(&temp);
 617                         }
 618
 619                         //append files
 620                         for(int i = 0 ; i < calculators.size(); i++) {
 621                                 string outputFile = outputDir + m->getRootName(m->getSimpleName(sharedfile)) + calculators[i]->getName() + ".amova";
 622
 623                                 for (int j = 0; j < processIDS.size(); j++) {
 624                                         m->appendFiles((outputFile + toString(processIDS[j])), outputFile);
 625                                         remove((outputFile + toString(processIDS[j])).c_str());
 626                                 }
 627                         }
 628                 }
 629 #else
 630                 driver(0, namesOfGroupCombos.size(), thisLookUp, "");
 631 #endif
 632
 633                 return 0;
 634         }
 635         catch(exception& e) {
 636                 m->errorOut(e, "AmovaCommand", "process");
 637                 exit(1);
 638         }
 639 }
 640 //**********************************************************************************************************************
 641
 642 int AmovaCommand::driver(int start, int num, vector<SharedRAbundVector*> thisLookup, string pidValue) {
 643         try {
 644                 vector<SharedRAbundVector*> subset;
 645                 EstOutput data;
 646
 647                 //for each calculator
 648                 for(int i = 0 ; i < calculators.size(); i++) {
 649
 650                         //create a new filename
 651                         ofstream out;
 652                         string outputFile = outputDir + m->getRootName(m->getSimpleName(sharedfile)) + calculators[i]->getName() + ".amova" + pidValue;
 653                         m->openOutputFileAppend(outputFile, out);
 654                         out.setf(ios::fixed, ios::floatfield); out.setf(ios::showpoint);
 655
 656                         //for each combo
 657                         for (int c = start; c < (start+num); c++) {
 658
 659                                 if (m->control_pressed) { out.close(); return 0; }
 660
 661                                 //get set names
 662                                 vector<string> setNames;
 663                                 for (int j = 0; j < namesOfGroupCombos[c].size(); j++) { setNames.push_back(namesOfGroupCombos[c][j]); }
 664
 665                                 vector<SharedRAbundVector*> thisCombosLookup;
 666                                 vector<string> thisCombosLookupSets; //what set each sharedRabund is from to be used when calculating SSWithin
 667                                 for (int k = 0; k < thisLookup.size(); k++) {
 668                                         string thisGroup = thisLookup[k]->getGroup();
 669
 670                                         //is this group for a set we want to compare??
 671                                         if (m->inUsersGroups(designMap->getGroup(thisGroup), setNames)) {
 672                                                 thisCombosLookup.push_back(thisLookup[k]);
 673                                                 thisCombosLookupSets.push_back(designMap->getGroup(thisGroup));
 674                                         }
 675
 676                                 }
 677
 678                                 int numGroups = thisCombosLookup.size();
 679
 680                                 //calc the distance matrix
 681                                 matrix.clear();
 682                                 matrix.resize(numGroups);
 683                                 for (int k = 0; k < matrix.size(); k++) {       for (int j = 0; j < matrix.size(); j++) {       matrix[k].push_back(1.0);       }       }
 684
 685                                 if (thisCombosLookup.size() == 0)  {
 686                                         m->mothurOut("[ERROR]: Missing shared info for sets. Skipping comparison."); m->mothurOutEndLine();
 687                                 }else{
 688
 689                                         out << thisLookup[0]->getLabel() << '\t';
 690                                         if (setNames.size() == 2) { out << setNames[0] << "-" << setNames[1] << '\t'; }
 691                                         else { out << "all" << '\t'; }
 692
 693                                         for (int k = 0; k < thisCombosLookup.size(); k++) {
 694                                                 for (int l = k; l < thisCombosLookup.size(); l++) {
 695
 696                                                         if (m->control_pressed) { out.close(); return 0; }
 697
 698                                                         if (k != l) { //we dont need to similiarity of a groups to itself
 699                                                                 //get estimated similarity between 2 groups
 700                                                                 subset.clear(); //clear out old pair of sharedrabunds
 701                                                                 //add new pair of sharedrabunds
 702                                                                 subset.push_back(thisCombosLookup[k]); subset.push_back(thisCombosLookup[l]);
 703
 704                                                                 //if this calc needs all groups to calculate the pair load all groups
 705                                                                 if (calculators[i]->getNeedsAll()) {
 706                                                                         //load subset with rest of lookup for those calcs that need everyone to calc for a pair
 707                                                                         for (int w = 0; w < thisCombosLookup.size(); w++) {
 708                                                                                 if ((w != k) && (w != l)) { subset.push_back(thisCombosLookup[w]); }
 709                                                                         }
 710                                                                 }
 711
 712                                                                 data = calculators[i]->getValues(subset); //saves the calculator outputs
 713
 714                                                                 //save values in similarity matrix
 715                                                                 matrix[k][l] = 1.0 - data[0];
 716                                                                 matrix[l][k] = 1.0 - data[0];
 717                                                         }
 718                                                 }
 719                                         }
 720
 721                                         //calc amova
 722                                         calcAmova(out, setNames.size(), thisCombosLookupSets);
 723                                 }
 724                         }
 725
 726                         out.close();
 727                 }
 728
 729                 return 0;
 730
 731         }
 732         catch(exception& e) {
 733                 m->errorOut(e, "AmovaCommand", "driver");
 734                 exit(1);
 735         }
 736 }
 737 //**********************************************************************************************************************
 738
 739 int AmovaCommand::driver(int start, int num, vector<string> names, string pidValue, vector< vector<double> >& completeMatrix) {
 740         try {
 741
 742                 //create a new filename
 743                 ofstream out;
 744                 string outputFile = outputDir + m->getRootName(m->getSimpleName(phylipfile)) + "amova" + pidValue;
 745                 m->openOutputFileAppend(outputFile, out);
 746                 out.setf(ios::fixed, ios::floatfield); out.setf(ios::showpoint);
 747
 748                 //for each combo
 749                 for (int c = start; c < (start+num); c++) {
 750
 751                         if (m->control_pressed) { out.close(); return 0; }
 752
 753                         //get set names
 754                         vector<string> setNames;
 755                         for (int j = 0; j < namesOfGroupCombos[c].size(); j++) { setNames.push_back(namesOfGroupCombos[c][j]); }
 756
 757                         vector<string> thisCombosSets; //what set each line in the distance matrix is from to be used when calculating SSWithin
 758                         set<int> indexes;
 759                         for (int k = 0; k < names.size(); k++) {
 760                                 //is this group for a set we want to compare??
 761                                 if (m->inUsersGroups(designMap->getGroup(names[k]), setNames)) {
 762                                         thisCombosSets.push_back(designMap->getGroup(names[k]));
 763                                         indexes.insert(k); //save indexes of valid rows in matrix for submatrix
 764                                 }
 765                         }
 766
 767                         int numGroups = thisCombosSets.size();
 768
 769                         //calc the distance matrix
 770                         matrix.clear();
 771                         matrix.resize(numGroups);
 772
 773                         for (int k = 0; k < matrix.size(); k++) {       for (int j = 0; j < matrix.size(); j++) {       matrix[k].push_back(1.0);       }       }
 774
 775                         if (thisCombosSets.size() == 0)  {
 776                                 m->mothurOut("[ERROR]: Missing distance info for sets. Skipping comparison."); m->mothurOutEndLine();
 777                         }else{
 778
 779                                 if (setNames.size() == 2) { out << setNames[0] << "-" << setNames[1] << '\t'; }
 780                                 else { out << "all" << '\t'; }
 781
 782                                 //fill submatrix
 783                                 int rowCount = 0;
 784                                 for (int j = 0; j < completeMatrix.size(); j++) {
 785
 786                                         if (indexes.count(j) != 0) { //we want at least part of this row
 787                                                 int columnCount = 0;
 788                                                 for (int k = 0; k < completeMatrix[j].size(); k++) {
 789
 790                                                         if (indexes.count(k) != 0) { //we want this distance
 791                                                                 matrix[rowCount][columnCount] = completeMatrix[j][k];
 792                                                                 matrix[columnCount][rowCount] = completeMatrix[j][k];
 793                                                                 columnCount++;
 794                                                         }
 795                                                 }
 796                                                 rowCount++;
 797                                         }
 798                                 }
 799
 800                                 //calc amova
 801                                 calcAmova(out, setNames.size(), thisCombosSets);
 802                         }
 803                 }
 804
 805                 out.close();
 806
 807
 808                 return 0;
 809
 810         }
 811         catch(exception& e) {
 812                 m->errorOut(e, "AmovaCommand", "driver");
 813                 exit(1);
 814         }
 815 }
 816 //**********************************************************************************************************************
 817 int AmovaCommand::calcAmova(ofstream& out, int numTreatments, vector<string> thisCombosLookupSets) {
 818         try {
 819
 820                 double SSWithin, SSTotal, SSAmoung, MSWithin, MSAmoung, FStatistic, pValue;
 821
 822                 SSWithin = calcWithin(matrix, numTreatments, thisCombosLookupSets);
 823                 SSTotal = calcTotal(numTreatments);
 824
 825                 int count = 0;
 826                 for (int i = 0; i < iters; i++) {
 827                         if (m->control_pressed) { break; }
 828
 829                         //randomly shuffle names to randomize the matrix
 830                         vector<string> copyNames = thisCombosLookupSets;
 831                         random_shuffle(copyNames.begin(), copyNames.end());
 832
 833                         double randomSSWithin = calcWithin(matrix, numTreatments, copyNames);
 834
 835                         if (randomSSWithin <= SSWithin) { count++; }
 836                 }
 837
 838                 SSAmoung = SSTotal - SSWithin;
 839                 MSWithin = SSWithin / (double) (matrix.size() - numTreatments);
 840                 MSAmoung = SSAmoung / (double) (numTreatments - 1);
 841                 FStatistic = MSAmoung / MSWithin;
 842
 843                 pValue = count / (float) iters;
 844
 845                 out << MSWithin << '\t' << MSAmoung << '\t' << FStatistic << '\t' << pValue << endl;
 846                 //cout << "ssAmong = " << SSAmoung << '\t' << "sswithin = " << SSWithin << '\t' << "ssTotal = " << SSTotal << '\t' << "pvalue = " << pValue << endl;
 847                 return 0;
 848         }
 849         catch(exception& e) {
 850                 m->errorOut(e, "AmovaCommand", "calcAmova");
 851                 exit(1);
 852         }
 853 }
 854 //**********************************************************************************************************************
 855 double AmovaCommand::calcWithin(vector< vector<double> >& thisMatrix, int numTreatments, vector<string> thisCombosLookupSets) {
 856         try {
 857                 double within = 0.0;
 858
 859                 //traverse lower triangle
 860                 for (int k = 0; k < thisMatrix.size(); k++) {
 861                         for (int l = k; l < thisMatrix[k].size(); l++) {
 862
 863                                 //if you are from the same treatment then eij is 1 so add, else eij = 0
 864                                 if (thisCombosLookupSets[k] == thisCombosLookupSets[l]) {
 865                                         within += (thisMatrix[k][l] * thisMatrix[k][l]); //dij^2
 866                                 }
 867                         }
 868                 }
 869
 870                 //1 / (numSamples / numTreatments)
 871                 within *= (1.0 / (float) (thisMatrix.size() / (float) numTreatments));
 872
 873                 return within;
 874         }
 875         catch(exception& e) {
 876                 m->errorOut(e, "AmovaCommand", "calcWithin");
 877                 exit(1);
 878         }
 879 }
 880 //**********************************************************************************************************************
 881 double AmovaCommand::calcTotal(int numTreatments) {
 882         try {
 883                 double total = 0.0;
 884
 885                 //traverse lower triangle
 886                 for (int k = 0; k < matrix.size(); k++) {
 887                         for (int l = k; l < matrix[k].size(); l++) {
 888                                 total += (matrix[k][l] * matrix[k][l]); //dij^2
 889                         }
 890                 }
 891
 892                 //1 / numSamples
 893                 total *= (1.0 / (float) matrix.size());
 894
 895                 return total;
 896         }
 897         catch(exception& e) {
 898                 m->errorOut(e, "AmovaCommand", "calcTotal");
 899                 exit(1);
 900         }
 901 }
 902 //**********************************************************************************************************************
 903