indicatorcommand.cpp

   1 /*
   2  *  indicatorcommand.cpp
   3  *  Mothur
   4  *
   5  *  Created by westcott on 11/12/10.
   6  *  Copyright 2010 Schloss Lab. All rights reserved.
   7  *
   8  */
   9
  10 #include "indicatorcommand.h"
  11 #include "sharedutilities.h"
  12
  13
  14 //**********************************************************************************************************************
  15 vector<string> IndicatorCommand::setParameters(){
  16         try {
  17                 CommandParameter piters("iters", "Number", "", "1000", "", "", "",false,false); parameters.push_back(piters);
  18                 CommandParameter pdesign("design", "InputTypes", "", "", "TreeDesign", "TreeDesign", "none",false,false); parameters.push_back(pdesign);
  19                 CommandParameter pshared("shared", "InputTypes", "", "", "SharedRel", "SharedRel", "none",false,false); parameters.push_back(pshared);
  20                 CommandParameter prelabund("relabund", "InputTypes", "", "", "SharedRel", "SharedRel", "none",false,false); parameters.push_back(prelabund);
  21                 CommandParameter pgroups("groups", "String", "", "", "", "", "",false,false); parameters.push_back(pgroups);
  22                 CommandParameter plabel("label", "String", "", "", "", "", "",false,false); parameters.push_back(plabel);
  23                 CommandParameter ptree("tree", "InputTypes", "", "", "TreeDesign", "TreeDesign", "none",false,false); parameters.push_back(ptree);
  24                 CommandParameter pinputdir("inputdir", "String", "", "", "", "", "",false,false); parameters.push_back(pinputdir);
  25                 CommandParameter poutputdir("outputdir", "String", "", "", "", "", "",false,false); parameters.push_back(poutputdir);
  26                 CommandParameter pprocessors("processors", "Number", "", "1", "", "", "",false,false); parameters.push_back(pprocessors);
  27
  28                 vector<string> myArray;
  29                 for (int i = 0; i < parameters.size(); i++) {   myArray.push_back(parameters[i].name);          }
  30                 return myArray;
  31         }
  32         catch(exception& e) {
  33                 m->errorOut(e, "IndicatorCommand", "setParameters");
  34                 exit(1);
  35         }
  36 }
  37 //**********************************************************************************************************************
  38 string IndicatorCommand::getHelpString(){
  39         try {
  40                 string helpString = "";
  41                 helpString += "The indicator command reads a shared or relabund file and a tree or design file, and outputs a .indicator.tre and .indicator.summary file. \n";
  42                 helpString += "The new tree contains labels at each internal node.  The label is the node number so you can relate the tree to the summary file.\n";
  43                 helpString += "The summary file lists the indicator value for each OTU for each node.\n";
  44                 helpString += "The indicator command parameters are tree, groups, shared, relabund, design and label. The tree parameter is required as well as either shared or relabund.\n";
  45                 helpString += "The design parameter allows you to relate the tree to the shared or relabund file, if your tree contains the grouping names, or if no tree is provided to group your groups into groupings.\n";
  46                 helpString += "The groups parameter allows you to specify which of the groups in your shared or relabund you would like analyzed, or if you provide a design file the groups in your design file.  The groups may be entered separated by dashes.\n";
  47                 helpString += "The label parameter indicates at what distance your tree relates to the shared or relabund.\n";
  48                 helpString += "The processors parameter allows you to specify how many processors you would like to use.  The default is 1. \n";
  49                 helpString += "The iters parameter allows you to set number of randomization for the P value.  The default is 1000.";
  50                 helpString += "The indicator command should be used in the following format: indicator(tree=test.tre, shared=test.shared, label=0.03)\n";
  51                 helpString += "Note: No spaces between parameter labels (i.e. tree), '=' and parameters (i.e.yourTreefile).\n";
  52                 return helpString;
  53         }
  54         catch(exception& e) {
  55                 m->errorOut(e, "IndicatorCommand", "getHelpString");
  56                 exit(1);
  57         }
  58 }
  59
  60 //**********************************************************************************************************************
  61 IndicatorCommand::IndicatorCommand(){
  62         try {
  63                 abort = true; calledHelp = true;
  64                 setParameters();
  65                 vector<string> tempOutNames;
  66                 outputTypes["tree"] = tempOutNames;
  67                 outputTypes["summary"] = tempOutNames;
  68         }
  69         catch(exception& e) {
  70                 m->errorOut(e, "IndicatorCommand", "IndicatorCommand");
  71                 exit(1);
  72         }
  73 }
  74 //**********************************************************************************************************************
  75 IndicatorCommand::IndicatorCommand(string option)  {
  76         try {
  77                 abort = false; calledHelp = false;
  78
  79                 //allow user to run help
  80                 if(option == "help") { help(); abort = true; calledHelp = true; }
  81                 else if(option == "citation") { citation(); abort = true; calledHelp = true;}
  82
  83                 else {
  84                         vector<string> myArray = setParameters();
  85
  86                         OptionParser parser(option);
  87                         map<string, string> parameters = parser.getParameters();
  88
  89                         ValidParameters validParameter;
  90                         map<string, string>::iterator it;
  91
  92                         //check to make sure all parameters are valid for command
  93                         for (it = parameters.begin(); it != parameters.end(); it++) {
  94                                 if (validParameter.isValidParameter(it->first, myArray, it->second) != true) {  abort = true;  }
  95                         }
  96
  97                         m->runParse = true;
  98                         m->Groups.clear();
  99                         m->namesOfGroups.clear();
 100                         m->Treenames.clear();
 101                         m->names.clear();
 102
 103                         vector<string> tempOutNames;
 104                         outputTypes["tree"] = tempOutNames;
 105                         outputTypes["summary"] = tempOutNames;
 106
 107                         //if the user changes the input directory command factory will send this info to us in the output parameter
 108                         string inputDir = validParameter.validFile(parameters, "inputdir", false);
 109                         if (inputDir == "not found"){   inputDir = "";          }
 110                         else {
 111                                 string path;
 112                                 it = parameters.find("tree");
 113                                 //user has given a template file
 114                                 if(it != parameters.end()){
 115                                         path = m->hasPath(it->second);
 116                                         //if the user has not given a path then, add inputdir. else leave path alone.
 117                                         if (path == "") {       parameters["tree"] = inputDir + it->second;             }
 118                                 }
 119
 120                                 it = parameters.find("shared");
 121                                 //user has given a template file
 122                                 if(it != parameters.end()){
 123                                         path = m->hasPath(it->second);
 124                                         //if the user has not given a path then, add inputdir. else leave path alone.
 125                                         if (path == "") {       parameters["shared"] = inputDir + it->second;           }
 126                                 }
 127
 128                                 it = parameters.find("relabund");
 129                                 //user has given a template file
 130                                 if(it != parameters.end()){
 131                                         path = m->hasPath(it->second);
 132                                         //if the user has not given a path then, add inputdir. else leave path alone.
 133                                         if (path == "") {       parameters["relabund"] = inputDir + it->second;         }
 134                                 }
 135
 136                                 it = parameters.find("design");
 137                                 //user has given a template file
 138                                 if(it != parameters.end()){
 139                                         path = m->hasPath(it->second);
 140                                         //if the user has not given a path then, add inputdir. else leave path alone.
 141                                         if (path == "") {       parameters["design"] = inputDir + it->second;           }
 142                                 }
 143                         }
 144
 145                         outputDir = validParameter.validFile(parameters, "outputdir", false);           if (outputDir == "not found"){  outputDir = ""; }
 146
 147                         //check for required parameters
 148                         treefile = validParameter.validFile(parameters, "tree", true);
 149                         if (treefile == "not open") { treefile = ""; abort = true; }
 150                         else if (treefile == "not found") { treefile = "";      }
 151                         else { m->setTreeFile(treefile); }
 152
 153                         sharedfile = validParameter.validFile(parameters, "shared", true);
 154                         if (sharedfile == "not open") { abort = true; }
 155                         else if (sharedfile == "not found") { sharedfile = ""; }
 156                         else { inputFileName = sharedfile; m->setSharedFile(sharedfile); }
 157
 158                         relabundfile = validParameter.validFile(parameters, "relabund", true);
 159                         if (relabundfile == "not open") { abort = true; }
 160                         else if (relabundfile == "not found") { relabundfile = ""; }
 161                         else { inputFileName = relabundfile; m->setRelAbundFile(relabundfile); }
 162
 163                         designfile = validParameter.validFile(parameters, "design", true);
 164                         if (designfile == "not open") { designfile = ""; abort = true; }
 165                         else if (designfile == "not found") { designfile = ""; }
 166                         else { m->setDesignFile(designfile); }
 167
 168                         groups = validParameter.validFile(parameters, "groups", false);
 169                         if (groups == "not found") { groups = "";  Groups.push_back("all"); }
 170                         else { m->splitAtDash(groups, Groups);  }
 171                         m->Groups = Groups;
 172
 173                         label = validParameter.validFile(parameters, "label", false);
 174                         if (label == "not found") { label = ""; m->mothurOut("You did not provide a label, I will use the first label in your inputfile."); m->mothurOutEndLine(); label=""; }
 175
 176                         string temp = validParameter.validFile(parameters, "iters", false);             if (temp == "not found") { temp = "1000"; }
 177                         convert(temp, iters);
 178
 179                         temp = validParameter.validFile(parameters, "processors", false);       if (temp == "not found"){       temp = m->getProcessors();      }
 180                         m->setProcessors(temp);
 181                         convert(temp, processors);
 182
 183                         if ((relabundfile == "") && (sharedfile == "")) {
 184                                 //is there are current file available for either of these?
 185                                 //give priority to shared, then relabund
 186                                 sharedfile = m->getSharedFile();
 187                                 if (sharedfile != "") {  inputFileName = sharedfile; m->mothurOut("Using " + sharedfile + " as input file for the shared parameter."); m->mothurOutEndLine(); }
 188                                 else {
 189                                         relabundfile = m->getRelAbundFile();
 190                                         if (relabundfile != "") {  inputFileName = relabundfile; m->mothurOut("Using " + relabundfile + " as input file for the relabund parameter."); m->mothurOutEndLine(); }
 191                                         else {
 192                                                 m->mothurOut("No valid current files. You must provide a shared or relabund."); m->mothurOutEndLine();
 193                                                 abort = true;
 194                                         }
 195                                 }
 196                         }
 197
 198
 199                         if ((designfile == "") && (treefile == "")) {
 200                                 treefile = m->getTreeFile();
 201                                 if (treefile != "") {  m->mothurOut("Using " + treefile + " as input file for the tree parameter."); m->mothurOutEndLine(); }
 202                                 else {
 203                                         designfile = m->getDesignFile();
 204                                         if (designfile != "") {  m->mothurOut("Using " + designfile + " as input file for the design parameter."); m->mothurOutEndLine(); }
 205                                         else {
 206                                                 m->mothurOut("[ERROR]: You must provide either a tree or design file."); m->mothurOutEndLine(); abort = true;
 207                                         }
 208                                 }
 209                         }
 210
 211                         if ((relabundfile != "") && (sharedfile != "")) { m->mothurOut("[ERROR]: You may not use both a shared and relabund file."); m->mothurOutEndLine(); abort = true;  }
 212
 213
 214                 }
 215         }
 216         catch(exception& e) {
 217                 m->errorOut(e, "IndicatorCommand", "IndicatorCommand");
 218                 exit(1);
 219         }
 220 }
 221 //**********************************************************************************************************************
 222
 223 int IndicatorCommand::execute(){
 224         try {
 225
 226                 if (abort == true) { if (calledHelp) { return 0; }  return 2;   }
 227
 228                 cout.setf(ios::fixed, ios::floatfield); cout.setf(ios::showpoint);
 229
 230                 int start = time(NULL);
 231
 232                 //read designfile if given and set up globaldatas groups for read of sharedfiles
 233                 if (designfile != "") {
 234                         designMap = new GroupMap(designfile);
 235                         designMap->readDesignMap();
 236
 237                         //fill Groups - checks for "all" and for any typo groups
 238                         SharedUtil* util = new SharedUtil();
 239                         util->setGroups(Groups, designMap->namesOfGroups);
 240                         delete util;
 241
 242                         //loop through the Groups and fill Globaldata's Groups with the design file info
 243                         m->Groups = designMap->getNamesSeqs(Groups);
 244                 }
 245
 246                 /***************************************************/
 247                 // use smart distancing to get right sharedRabund  //
 248                 /***************************************************/
 249                 if (sharedfile != "") {
 250                         getShared();
 251                         if (m->control_pressed) { if (designfile != "") { delete designMap; } for (int i = 0; i < lookup.size(); i++) {  delete lookup[i];  } return 0; }
 252                         if (lookup[0] == NULL) { m->mothurOut("[ERROR] reading shared file."); m->mothurOutEndLine(); return 0; }
 253                 }else {
 254                         getSharedFloat();
 255                         if (m->control_pressed) {  if (designfile != "") { delete designMap; } for (int i = 0; i < lookupFloat.size(); i++) {  delete lookupFloat[i];  } return 0; }
 256                         if (lookupFloat[0] == NULL) { m->mothurOut("[ERROR] reading relabund file."); m->mothurOutEndLine(); return 0; }
 257                 }
 258
 259                 //reset Globaldatas groups if needed
 260                 if (designfile != "") { m->Groups = Groups; }
 261
 262                 /***************************************************/
 263                 //    reading tree info                                                    //
 264                 /***************************************************/
 265                 if (treefile != "") {
 266                         string groupfile = "";
 267                         m->setTreeFile(treefile);
 268                         Tree* tree = new Tree(treefile); delete tree;  //extracts names from tree to make faked out groupmap
 269                         treeMap = new TreeMap();
 270                         bool mismatch = false;
 271
 272                         for (int i = 0; i < m->Treenames.size(); i++) {
 273                                 //sanity check - is this a group that is not in the sharedfile?
 274                                 if (designfile == "") {
 275                                         if (!(m->inUsersGroups(m->Treenames[i], m->namesOfGroups))) {
 276                                                 m->mothurOut("[ERROR]: " + m->Treenames[i] + " is not a group in your shared or relabund file."); m->mothurOutEndLine();
 277                                                 mismatch = true;
 278                                         }
 279                                         treeMap->addSeq(m->Treenames[i], "Group1");
 280                                 }else{
 281                                         vector<string> myGroups; myGroups.push_back(m->Treenames[i]);
 282                                         vector<string> myNames = designMap->getNamesSeqs(myGroups);
 283
 284                                         for(int k = 0; k < myNames.size(); k++) {
 285                                                 if (!(m->inUsersGroups(myNames[k], m->namesOfGroups))) {
 286                                                         m->mothurOut("[ERROR]: " + myNames[k] + " is not a group in your shared or relabund file."); m->mothurOutEndLine();
 287                                                         mismatch = true;
 288                                                 }
 289                                         }
 290                                         treeMap->addSeq(m->Treenames[i], "Group1");
 291                                 }
 292                         }
 293
 294                         if ((designfile != "") && (m->Treenames.size() != Groups.size())) { cout << Groups.size() << '\t' << m->Treenames.size() << endl; m->mothurOut("[ERROR]: You design file does not match your tree, aborting."); m->mothurOutEndLine(); mismatch = true; }
 295
 296                         if (mismatch) { //cleanup and exit
 297                                 if (designfile != "") { delete designMap; }
 298                                 if (sharedfile != "") {  for (int i = 0; i < lookup.size(); i++) {  delete lookup[i];  } }
 299                                 else { for (int i = 0; i < lookupFloat.size(); i++) {  delete lookupFloat[i];  } }
 300                                 delete treeMap;
 301                                 return 0;
 302                         }
 303
 304                         read = new ReadNewickTree(treefile);
 305                         int readOk = read->read(treeMap);
 306
 307                         if (readOk != 0) { m->mothurOut("Read Terminated."); m->mothurOutEndLine(); delete treeMap; delete read; return 0; }
 308
 309                         vector<Tree*> T = read->getTrees();
 310
 311                         delete read;
 312
 313                         if (m->control_pressed) {
 314                                 if (designfile != "") { delete designMap; }
 315                                 if (sharedfile != "") {  for (int i = 0; i < lookup.size(); i++) {  delete lookup[i];  } }
 316                                 else { for (int i = 0; i < lookupFloat.size(); i++) {  delete lookupFloat[i];  } }
 317                                 for (int i = 0; i < T.size(); i++) {  delete T[i];  }  delete treeMap; return 0;
 318                         }
 319
 320                         T[0]->assembleTree();
 321
 322                         /***************************************************/
 323                         //    create ouptut tree - respecting pickedGroups //
 324                         /***************************************************/
 325                         Tree* outputTree = new Tree(m->Groups.size(), treeMap);
 326
 327                         outputTree->getSubTree(T[0], m->Groups);
 328                         outputTree->assembleTree();
 329
 330                         //no longer need original tree, we have output tree to use and label
 331                         for (int i = 0; i < T.size(); i++) {  delete T[i];  }
 332
 333
 334                         if (m->control_pressed) {
 335                                 if (designfile != "") { delete designMap; }
 336                                 if (sharedfile != "") {  for (int i = 0; i < lookup.size(); i++) {  delete lookup[i];  } }
 337                                 else { for (int i = 0; i < lookupFloat.size(); i++) {  delete lookupFloat[i];  } }
 338                                 delete outputTree; delete treeMap;  return 0;
 339                         }
 340
 341                         /***************************************************/
 342                         //              get indicator species values                       //
 343                         /***************************************************/
 344                         GetIndicatorSpecies(outputTree);
 345                         delete outputTree; delete treeMap;
 346
 347                 }else { //run with design file only
 348                         //get indicator species
 349                         GetIndicatorSpecies();
 350                 }
 351
 352                 if (designfile != "") { delete designMap; }
 353                 if (sharedfile != "") {  for (int i = 0; i < lookup.size(); i++) {  delete lookup[i];  } }
 354                 else { for (int i = 0; i < lookupFloat.size(); i++) {  delete lookupFloat[i];  } }
 355
 356                 if (m->control_pressed) { for (int i = 0; i < outputNames.size(); i++) {        remove(outputNames[i].c_str()); } return 0; }
 357
 358                 //set tree file as new current treefile
 359                 if (treefile != "") {
 360                         string current = "";
 361                         itTypes = outputTypes.find("tree");
 362                         if (itTypes != outputTypes.end()) {
 363                                 if ((itTypes->second).size() != 0) { current = (itTypes->second)[0]; m->setTreeFile(current); }
 364                         }
 365                 }
 366
 367                 m->mothurOutEndLine(); m->mothurOutEndLine(); m->mothurOut("It took " + toString(time(NULL) - start) + " secs to find the indicator species."); m->mothurOutEndLine();
 368                 m->mothurOutEndLine();
 369                 m->mothurOut("Output File Names: "); m->mothurOutEndLine();
 370                 for (int i = 0; i < outputNames.size(); i++) {  m->mothurOut(outputNames[i]); m->mothurOutEndLine();    }
 371                 m->mothurOutEndLine();
 372
 373                 return 0;
 374         }
 375         catch(exception& e) {
 376                 m->errorOut(e, "IndicatorCommand", "execute");
 377                 exit(1);
 378         }
 379 }
 380 //**********************************************************************************************************************
 381 //divide shared or relabund file by groupings in the design file
 382 //report all otu values to file
 383 int IndicatorCommand::GetIndicatorSpecies(){
 384         try {
 385                 string thisOutputDir = outputDir;
 386                 if (outputDir == "") {  thisOutputDir += m->hasPath(inputFileName);  }
 387                 string outputFileName = thisOutputDir + m->getRootName(m->getSimpleName(inputFileName)) + "indicator.summary";
 388                 outputNames.push_back(outputFileName); outputTypes["summary"].push_back(outputFileName);
 389
 390                 ofstream out;
 391                 m->openOutputFile(outputFileName, out);
 392                 out.setf(ios::fixed, ios::floatfield); out.setf(ios::showpoint);
 393                 m->mothurOutEndLine(); m->mothurOut("Species\tIndicatorValue\tpValue\n");
 394
 395                 int numBins = 0;
 396                 if (sharedfile != "") { numBins = lookup[0]->getNumBins(); }
 397                 else { numBins = lookupFloat[0]->getNumBins(); }
 398
 399                 if (m->control_pressed) { out.close(); return 0; }
 400
 401                 /*****************************************************/
 402                 //create vectors containing rabund info                          //
 403                 /*****************************************************/
 404
 405                 vector<float> indicatorValues; //size of numBins
 406                 vector<float> pValues;
 407                 map< vector<int>, vector<int> > randomGroupingsMap; //maps location in groupings to location in groupings, ie, [0][0] -> [1][2]. This is so we don't have to actually move the sharedRabundVectors.
 408
 409                 if (sharedfile != "") {
 410                         vector< vector<SharedRAbundVector*> > groupings;
 411                         set<string> groupsAlreadyAdded;
 412                         vector<SharedRAbundVector*> subset;
 413
 414                         //for each grouping
 415                         for (int i = 0; i < designMap->namesOfGroups.size(); i++) {
 416
 417                                 for (int k = 0; k < lookup.size(); k++) {
 418                                         //are you from this grouping?
 419                                         if (designMap->getGroup(lookup[k]->getGroup()) == designMap->namesOfGroups[i]) {
 420                                                 subset.push_back(lookup[k]);
 421                                                 groupsAlreadyAdded.insert(lookup[k]->getGroup());
 422                                         }
 423                                 }
 424                                 if (subset.size() != 0) { groupings.push_back(subset); }
 425                                 subset.clear();
 426                         }
 427
 428                         if (groupsAlreadyAdded.size() != lookup.size()) {  m->mothurOut("[ERROR]: could not make proper groupings."); m->mothurOutEndLine(); }
 429
 430                         indicatorValues = getValues(groupings, randomGroupingsMap);
 431
 432                         pValues = getPValues(groupings, randomGroupingsMap, lookup.size(), indicatorValues);
 433                 }else {
 434                         vector< vector<SharedRAbundFloatVector*> > groupings;
 435                         set<string> groupsAlreadyAdded;
 436                         vector<SharedRAbundFloatVector*> subset;
 437
 438                         //for each grouping
 439                         for (int i = 0; i < designMap->namesOfGroups.size(); i++) {
 440                                 for (int k = 0; k < lookupFloat.size(); k++) {
 441                                         //are you from this grouping?
 442                                         if (designMap->getGroup(lookupFloat[k]->getGroup()) == designMap->namesOfGroups[i]) {
 443                                                 subset.push_back(lookupFloat[k]);
 444                                                 groupsAlreadyAdded.insert(lookupFloat[k]->getGroup());
 445                                         }
 446                                 }
 447                                 if (subset.size() != 0) { groupings.push_back(subset); }
 448                                 subset.clear();
 449                         }
 450
 451                         if (groupsAlreadyAdded.size() != lookupFloat.size()) {  m->mothurOut("[ERROR]: could not make proper groupings."); m->mothurOutEndLine(); }
 452
 453                         indicatorValues = getValues(groupings, randomGroupingsMap);
 454
 455                         pValues = getPValues(groupings, randomGroupingsMap, lookupFloat.size(), indicatorValues);
 456                 }
 457
 458                 if (m->control_pressed) { out.close(); return 0; }
 459
 460
 461                 /******************************************************/
 462                 //output indicator values to table form               //
 463                 /*****************************************************/
 464                 out << "OTU\tIndicator_Value\tpValue" << endl;
 465                 for (int j = 0; j < indicatorValues.size(); j++) {
 466
 467                         if (m->control_pressed) { out.close(); return 0; }
 468
 469                         out << (j+1) << '\t' << indicatorValues[j] << '\t';
 470
 471                         if (pValues[j] > (1/(float)iters)) { out << pValues[j] << endl; }
 472                         else { out << "<" << (1/(float)iters) << endl; }
 473
 474                         if (pValues[j] <= 0.05) {
 475                                 cout << "OTU" << j+1 << '\t' << indicatorValues[j]  << '\t';
 476                                 string pValueString = "<" + toString((1/(float)iters));
 477                                 if (pValues[j] > (1/(float)iters)) { pValueString = toString(pValues[j]); cout << pValues[j];}
 478                                 else { cout << "<" << (1/(float)iters); }
 479                                 m->mothurOutJustToLog("OTU" + toString(j+1) + "\t" + toString(indicatorValues[j]) + "\t" + pValueString);
 480                                 m->mothurOutEndLine();
 481                         }
 482                 }
 483
 484                 out.close();
 485
 486                 return 0;
 487         }
 488         catch(exception& e) {
 489                 m->errorOut(e, "IndicatorCommand", "GetIndicatorSpecies");
 490                 exit(1);
 491         }
 492 }
 493 //**********************************************************************************************************************
 494 //traverse tree finding indicator species values for each otu at each node
 495 //label node with otu number that has highest indicator value
 496 //report all otu values to file
 497 int IndicatorCommand::GetIndicatorSpecies(Tree*& T){
 498         try {
 499
 500                 string thisOutputDir = outputDir;
 501                 if (outputDir == "") {  thisOutputDir += m->hasPath(inputFileName);  }
 502                 string outputFileName = thisOutputDir + m->getRootName(m->getSimpleName(inputFileName)) + "indicator.summary";
 503                 outputNames.push_back(outputFileName); outputTypes["summary"].push_back(outputFileName);
 504
 505                 ofstream out;
 506                 m->openOutputFile(outputFileName, out);
 507                 out.setf(ios::fixed, ios::floatfield); out.setf(ios::showpoint);
 508
 509                 int numBins = 0;
 510                 if (sharedfile != "") { numBins = lookup[0]->getNumBins(); }
 511                 else { numBins = lookupFloat[0]->getNumBins(); }
 512
 513                 //print headings
 514                 out << "TreeNode\t";
 515                 for (int i = 0; i < numBins; i++) { out << "OTU" << (i+1) << "_IndValue" << '\t' << "pValue" << '\t'; }
 516                 out << endl;
 517
 518                 m->mothurOutEndLine(); m->mothurOut("Node\tSpecies\tIndicatorValue\tpValue\n");
 519
 520                 string treeOutputDir = outputDir;
 521                 if (outputDir == "") {  treeOutputDir += m->hasPath(treefile);  }
 522                 string outputTreeFileName = treeOutputDir + m->getRootName(m->getSimpleName(treefile)) + "indicator.tre";
 523
 524
 525                 //create a map from tree node index to names of descendants, save time later to know which sharedRabund you need
 526                 map<int, set<string> > nodeToDescendants;
 527                 map<int, set<int> > descendantNodes;
 528                 for (int i = 0; i < T->getNumNodes(); i++) {
 529                         if (m->control_pressed) { return 0; }
 530
 531                         nodeToDescendants[i] = getDescendantList(T, i, nodeToDescendants, descendantNodes);
 532                 }
 533
 534                 //you need the distances to leaf to decide grouping below
 535                 //this will also set branch lengths if the tree does not include them
 536                 map<int, float> distToRoot = getDistToRoot(T);
 537
 538                 //for each node
 539                 for (int i = T->getNumLeaves(); i < T->getNumNodes(); i++) {
 540                         //cout << endl << i+1 << endl;
 541                         if (m->control_pressed) { out.close(); return 0; }
 542
 543                         /*****************************************************/
 544                         //create vectors containing rabund info                          //
 545                         /*****************************************************/
 546
 547                         vector<float> indicatorValues; //size of numBins
 548                         vector<float> pValues;
 549                         map< vector<int>, vector<int> > randomGroupingsMap; //maps location in groupings to location in groupings, ie, [0][0] -> [1][2]. This is so we don't have to actually move the sharedRabundVectors.
 550
 551                         if (sharedfile != "") {
 552                                 vector< vector<SharedRAbundVector*> > groupings;
 553
 554                                 //get nodes that will be a valid grouping
 555                                 //you are valid if you are not one of my descendants
 556                                 //AND your distToRoot is >= mine
 557                                 //AND you were not added as part of a larger grouping. Largest nodes are added first.
 558
 559                                 set<string> groupsAlreadyAdded;
 560                                 //create a grouping with my grouping
 561                                 vector<SharedRAbundVector*> subset;
 562                                 int count = 0;
 563                                 int doneCount = nodeToDescendants[i].size();
 564                                 for (int k = 0; k < lookup.size(); k++) {
 565                                         //is this descendant of i
 566                                         if ((nodeToDescendants[i].count(lookup[k]->getGroup()) != 0)) {
 567                                                 subset.push_back(lookup[k]);
 568                                                 groupsAlreadyAdded.insert(lookup[k]->getGroup());
 569                                                 count++;
 570                                         }
 571                                         if (count == doneCount) { break; } //quit once you get the rabunds for this grouping
 572                                 }
 573                                 if (subset.size() != 0) { groupings.push_back(subset); }
 574
 575
 576                                 for (int j = (T->getNumNodes()-1); j >= 0; j--) {
 577
 578
 579                                         if ((descendantNodes[i].count(j) == 0) && (distToRoot[j] >= distToRoot[i])) {
 580                                                 vector<SharedRAbundVector*> subset;
 581                                                 int count = 0;
 582                                                 int doneCount = nodeToDescendants[j].size();
 583                                                 for (int k = 0; k < lookup.size(); k++) {
 584                                                         //is this descendant of j, and we didn't already add this as part of a larger grouping
 585                                                         if ((nodeToDescendants[j].count(lookup[k]->getGroup()) != 0) && (groupsAlreadyAdded.count(lookup[k]->getGroup()) == 0)) {
 586                                                                 subset.push_back(lookup[k]);
 587                                                                 groupsAlreadyAdded.insert(lookup[k]->getGroup());
 588                                                                 count++;
 589                                                         }
 590                                                         if (count == doneCount) { break; } //quit once you get the rabunds for this grouping
 591                                                 }
 592
 593                                                 //if subset.size == 0 then the node was added as part of a larger grouping
 594                                                 if (subset.size() != 0) { groupings.push_back(subset); }
 595                                         }
 596                                 }
 597
 598                                 if (groupsAlreadyAdded.size() != lookup.size()) {  m->mothurOut("[ERROR]: could not make proper groupings."); m->mothurOutEndLine(); }
 599
 600                                 indicatorValues = getValues(groupings, randomGroupingsMap);
 601
 602                                 pValues = getPValues(groupings, randomGroupingsMap, lookup.size(), indicatorValues);
 603                         }else {
 604                                 vector< vector<SharedRAbundFloatVector*> > groupings;
 605
 606                                 //get nodes that will be a valid grouping
 607                                 //you are valid if you are not one of my descendants
 608                                 //AND your distToRoot is >= mine
 609                                 //AND you were not added as part of a larger grouping. Largest nodes are added first.
 610
 611                                 set<string> groupsAlreadyAdded;
 612                                 //create a grouping with my grouping
 613                                 vector<SharedRAbundFloatVector*> subset;
 614                                 int count = 0;
 615                                 int doneCount = nodeToDescendants[i].size();
 616                                 for (int k = 0; k < lookupFloat.size(); k++) {
 617                                         //is this descendant of i
 618                                         if ((nodeToDescendants[i].count(lookupFloat[k]->getGroup()) != 0)) {
 619                                                 subset.push_back(lookupFloat[k]);
 620                                                 groupsAlreadyAdded.insert(lookupFloat[k]->getGroup());
 621                                                 count++;
 622                                         }
 623                                         if (count == doneCount) { break; } //quit once you get the rabunds for this grouping
 624                                 }
 625                                 if (subset.size() != 0) { groupings.push_back(subset); }
 626
 627                                 for (int j = (T->getNumNodes()-1); j >= 0; j--) {
 628                                         if ((descendantNodes[i].count(j) == 0) && (distToRoot[j] >= distToRoot[i])) {
 629                                                 vector<SharedRAbundFloatVector*> subset;
 630                                                 int count = 0;
 631                                                 int doneCount = nodeToDescendants[j].size();
 632                                                 for (int k = 0; k < lookupFloat.size(); k++) {
 633                                                         //is this descendant of j, and we didn't already add this as part of a larger grouping
 634                                                         if ((nodeToDescendants[j].count(lookupFloat[k]->getGroup()) != 0) && (groupsAlreadyAdded.count(lookupFloat[k]->getGroup()) == 0)) {
 635                                                                 subset.push_back(lookupFloat[k]);
 636                                                                 groupsAlreadyAdded.insert(lookupFloat[k]->getGroup());
 637                                                                 count++;
 638                                                         }
 639                                                         if (count == doneCount) { break; } //quit once you get the rabunds for this grouping
 640                                                 }
 641
 642                                                 //if subset.size == 0 then the node was added as part of a larger grouping
 643                                                 if (subset.size() != 0) { groupings.push_back(subset); }
 644                                         }
 645                                 }
 646
 647                                 if (groupsAlreadyAdded.size() != lookupFloat.size()) { m->mothurOut("[ERROR]: could not make proper groupings."); m->mothurOutEndLine(); }
 648
 649                                 indicatorValues = getValues(groupings, randomGroupingsMap);
 650
 651                                 pValues = getPValues(groupings, randomGroupingsMap, lookupFloat.size(), indicatorValues);
 652                         }
 653
 654                         if (m->control_pressed) { out.close(); return 0; }
 655
 656
 657                         /******************************************************/
 658                         //output indicator values to table form + label tree  //
 659                         /*****************************************************/
 660                         out << (i+1) << '\t';
 661                         for (int j = 0; j < indicatorValues.size(); j++) {
 662
 663                                 if (m->control_pressed) { out.close(); return 0; }
 664
 665                                 if (pValues[j] < (1/(float)iters)) {
 666                                         out << indicatorValues[j] << '\t' << '<' << (1/(float)iters) << '\t';
 667                                 }else {
 668                                         out << indicatorValues[j] << '\t' << pValues[j] << '\t';
 669                                 }
 670
 671                                 if (pValues[j] <= 0.05) {
 672                                         cout << i+1 << "\tOTU" << j+1 << '\t' << indicatorValues[j]  << '\t';
 673                                         string pValueString = "<" + toString((1/(float)iters));
 674                                         if (pValues[j] > (1/(float)iters)) { pValueString = toString(pValues[j]); cout << pValues[j];}
 675                                         else { cout << "<" << (1/(float)iters); }
 676                                         m->mothurOutJustToLog(toString(i) + "\tOTU" + toString(j+1) + "\t" + toString(indicatorValues[j]) + "\t" + pValueString);
 677                                         m->mothurOutEndLine();
 678                                 }
 679                         }
 680                         out << endl;
 681
 682                         T->tree[i].setLabel((i+1));
 683                 }
 684                 out.close();
 685
 686                 ofstream outTree;
 687                 m->openOutputFile(outputTreeFileName, outTree);
 688                 outputNames.push_back(outputTreeFileName); outputTypes["tree"].push_back(outputTreeFileName);
 689
 690                 T->print(outTree, "both");
 691                 outTree.close();
 692
 693                 return 0;
 694         }
 695         catch(exception& e) {
 696                 m->errorOut(e, "IndicatorCommand", "GetIndicatorSpecies");
 697                 exit(1);
 698         }
 699 }
 700 //**********************************************************************************************************************
 701 vector<float> IndicatorCommand::getValues(vector< vector<SharedRAbundFloatVector*> >& groupings, map< vector<int>, vector<int> > groupingsMap){
 702         try {
 703                 vector<float> values;
 704                 map< vector<int>, vector<int> >::iterator it;
 705
 706                 //for each otu
 707                 for (int i = 0; i < groupings[0][0]->getNumBins(); i++) {
 708
 709                         if (m->control_pressed) { return values; }
 710
 711                         vector<float> terms;
 712                         float AijDenominator = 0.0;
 713                         vector<float> Bij;
 714
 715                         //get overall abundance of each grouping
 716                         for (int j = 0; j < groupings.size(); j++) {
 717
 718                                 float totalAbund = 0;
 719                                 int numNotZero = 0;
 720                                 for (int k = 0; k < groupings[j].size(); k++) {
 721                                         vector<int> temp; temp.push_back(j); temp.push_back(k);
 722                                         it = groupingsMap.find(temp);
 723
 724                                         if (it == groupingsMap.end()) { //this one didnt get moved
 725                                                 totalAbund += groupings[j][k]->getAbundance(i);
 726                                                 if (groupings[j][k]->getAbundance(i) != 0.0) { numNotZero++; }
 727                                         }else {
 728                                                 totalAbund += groupings[(it->second)[0]][(it->second)[1]]->getAbundance(i);
 729                                                 if (groupings[(it->second)[0]][(it->second)[1]]->getAbundance(i) != 0.0) { numNotZero++; }
 730                                         }
 731
 732                                 }
 733
 734                                 //mean abundance
 735                                 float Aij = (totalAbund / (float) groupings[j].size());
 736                                 terms.push_back(Aij);
 737
 738                                 //percentage of sites represented
 739                                 Bij.push_back(numNotZero / (float) groupings[j].size());
 740
 741                                 AijDenominator += Aij;
 742                         }
 743
 744                         float maxIndVal = 0.0;
 745                         for (int j = 0; j < terms.size(); j++) {
 746                                 float thisAij = (terms[j] / AijDenominator); //relative abundance
 747                                 float thisValue = thisAij * Bij[j] * 100.0;
 748
 749                                 //save largest
 750                                 if (thisValue > maxIndVal) { maxIndVal = thisValue; }
 751                         }
 752
 753                         values.push_back(maxIndVal);
 754                 }
 755
 756                 return values;
 757         }
 758         catch(exception& e) {
 759                 m->errorOut(e, "IndicatorCommand", "getValues");
 760                 exit(1);
 761         }
 762 }
 763 //**********************************************************************************************************************
 764 //same as above, just data type difference
 765 vector<float> IndicatorCommand::getValues(vector< vector<SharedRAbundVector*> >& groupings, map< vector<int>, vector<int> > groupingsMap){
 766         try {
 767                 vector<float> values;
 768
 769         /*for (int j = 0; j < groupings.size(); j++) {
 770                 cout << "grouping " << j << endl;
 771                 for (int k = 0; k < groupings[j].size(); k++) {
 772                         cout << groupings[j][k]->getGroup() << endl;
 773                 }
 774         }*/
 775                 map< vector<int>, vector<int> >::iterator it;
 776
 777                 //for each otu
 778                 for (int i = 0; i < groupings[0][0]->getNumBins(); i++) {
 779                         vector<float> terms;
 780                         float AijDenominator = 0.0;
 781                         vector<float> Bij;
 782                         //get overall abundance of each grouping
 783                         for (int j = 0; j < groupings.size(); j++) {
 784
 785                                 int totalAbund = 0.0;
 786                                 int numNotZero = 0;
 787                                 for (int k = 0; k < groupings[j].size(); k++) {
 788                                         vector<int> temp; temp.push_back(j); temp.push_back(k);
 789                                         it = groupingsMap.find(temp);
 790
 791                                         if (it == groupingsMap.end()) { //this one didnt get moved
 792                                                 totalAbund += groupings[j][k]->getAbundance(i);
 793                                                 if (groupings[j][k]->getAbundance(i) != 0.0) { numNotZero++; }
 794                                         }else {
 795                                                 totalAbund += groupings[(it->second)[0]][(it->second)[1]]->getAbundance(i);
 796                                                 if (groupings[(it->second)[0]][(it->second)[1]]->getAbundance(i) != 0.0) { numNotZero++; }
 797                                         }
 798
 799                                 }
 800
 801                                 //mean abundance
 802                                 float Aij = (totalAbund / (float) groupings[j].size());
 803                                 terms.push_back(Aij);
 804
 805                                 //percentage of sites represented
 806                                 Bij.push_back(numNotZero / (float) groupings[j].size());
 807
 808                                 AijDenominator += Aij;
 809                         }
 810
 811                         float maxIndVal = 0.0;
 812                         for (int j = 0; j < terms.size(); j++) {
 813                                 float thisAij = (terms[j] / AijDenominator); //relative abundance
 814                                 float thisValue = thisAij * Bij[j] * 100.0;
 815
 816                                 //save largest
 817                                 if (thisValue > maxIndVal) { maxIndVal = thisValue; }
 818                         }
 819
 820                         values.push_back(maxIndVal);
 821                 }
 822
 823                 return values;
 824         }
 825         catch(exception& e) {
 826                 m->errorOut(e, "IndicatorCommand", "getValues");
 827                 exit(1);
 828         }
 829 }
 830 //**********************************************************************************************************************
 831 //you need the distances to root to decide groupings
 832 //this will also set branch lengths if the tree does not include them
 833 map<int, float> IndicatorCommand::getDistToRoot(Tree*& T){
 834         try {
 835                 map<int, float> dists;
 836
 837                 bool hasBranchLengths = false;
 838                 for (int i = 0; i < T->getNumNodes(); i++) {
 839                         if (T->tree[i].getBranchLength() > 0.0) {  hasBranchLengths = true; break; }
 840                 }
 841
 842                 //set branchlengths if needed
 843                 if (!hasBranchLengths) {
 844                         for (int i = 0; i < T->getNumNodes(); i++) {
 845
 846                                 int lc = T->tree[i].getLChild();
 847                                 int rc = T->tree[i].getRChild();
 848
 849                                 if (lc == -1) { // you are a leaf
 850                                         //if you are a leaf set you priliminary length to 1.0, this may adjust later
 851                                         T->tree[i].setBranchLength(1.0);
 852                                         dists[i] = 1.0;
 853                                 }else{ // you are an internal node
 854                                         //look at your children's length to leaf
 855                                         float ldist = dists[lc];
 856                                         float rdist = dists[rc];
 857
 858                                         float greater = ldist;
 859                                         if (rdist > greater) { greater = rdist; dists[i] = ldist + 1.0;}
 860                                         else { dists[i] = rdist + 1.0; }
 861
 862
 863                                         //branch length = difference + 1
 864                                         T->tree[lc].setBranchLength((abs(ldist-greater) + 1.0));
 865                                         T->tree[rc].setBranchLength((abs(rdist-greater) + 1.0));
 866                                 }
 867                         }
 868                 }
 869
 870                 dists.clear();
 871
 872                 for (int i = 0; i < T->getNumNodes(); i++) {
 873
 874                         double sum = 0.0;
 875                         int index = i;
 876
 877                         while(T->tree[index].getParent() != -1){
 878                                 if (T->tree[index].getBranchLength() != -1) {
 879                                         sum += abs(T->tree[index].getBranchLength());
 880                                 }
 881                                 index = T->tree[index].getParent();
 882                         }
 883
 884                         dists[i] = sum;
 885                 }
 886
 887                 return dists;
 888         }
 889         catch(exception& e) {
 890                 m->errorOut(e, "IndicatorCommand", "getLengthToLeaf");
 891                 exit(1);
 892         }
 893 }
 894 //**********************************************************************************************************************
 895 set<string> IndicatorCommand::getDescendantList(Tree*& T, int i, map<int, set<string> > descendants, map<int, set<int> >& nodes){
 896         try {
 897                 set<string> names;
 898
 899                 set<string>::iterator it;
 900
 901                 int lc = T->tree[i].getLChild();
 902                 int rc = T->tree[i].getRChild();
 903
 904                 if (lc == -1) { //you are a leaf your only descendant is yourself
 905                         set<int> temp; temp.insert(i);
 906                         nodes[i] = temp;
 907
 908                         if (designfile == "") {
 909                                 names.insert(T->tree[i].getName());
 910                         }else {
 911                                 vector<string> myGroup; myGroup.push_back(T->tree[i].getName());
 912                                 vector<string> myReps = designMap->getNamesSeqs(myGroup);
 913                                 for (int k = 0; k < myReps.size(); k++) {
 914                                         names.insert(myReps[k]);
 915                                 }
 916                         }
 917
 918                 }else{ //your descedants are the combination of your childrens descendants
 919                         names = descendants[lc];
 920                         nodes[i] = nodes[lc];
 921                         for (it = descendants[rc].begin(); it != descendants[rc].end(); it++) {
 922                                 names.insert(*it);
 923                         }
 924                         for (set<int>::iterator itNum = nodes[rc].begin(); itNum != nodes[rc].end(); itNum++) {
 925                                 nodes[i].insert(*itNum);
 926                         }
 927                         //you are your own descendant
 928                         nodes[i].insert(i);
 929                 }
 930
 931                 return names;
 932         }
 933         catch(exception& e) {
 934                 m->errorOut(e, "IndicatorCommand", "getDescendantList");
 935                 exit(1);
 936         }
 937 }
 938 //**********************************************************************************************************************
 939 int IndicatorCommand::getShared(){
 940         try {
 941                 InputData* input = new InputData(sharedfile, "sharedfile");
 942                 lookup = input->getSharedRAbundVectors();
 943                 string lastLabel = lookup[0]->getLabel();
 944
 945                 if (label == "") { label = lastLabel; delete input; return 0; }
 946
 947                 //if the users enters label "0.06" and there is no "0.06" in their file use the next lowest label.
 948                 set<string> labels; labels.insert(label);
 949                 set<string> processedLabels;
 950                 set<string> userLabels = labels;
 951
 952                 //as long as you are not at the end of the file or done wih the lines you want
 953                 while((lookup[0] != NULL) && (userLabels.size() != 0)) {
 954                         if (m->control_pressed) {  delete input; return 0;  }
 955
 956                         if(labels.count(lookup[0]->getLabel()) == 1){
 957                                 processedLabels.insert(lookup[0]->getLabel());
 958                                 userLabels.erase(lookup[0]->getLabel());
 959                                 break;
 960                         }
 961
 962                         if ((m->anyLabelsToProcess(lookup[0]->getLabel(), userLabels, "") == true) && (processedLabels.count(lastLabel) != 1)) {
 963                                 string saveLabel = lookup[0]->getLabel();
 964
 965                                 for (int i = 0; i < lookup.size(); i++) {  delete lookup[i];  }
 966                                 lookup = input->getSharedRAbundVectors(lastLabel);
 967
 968                                 processedLabels.insert(lookup[0]->getLabel());
 969                                 userLabels.erase(lookup[0]->getLabel());
 970
 971                                 //restore real lastlabel to save below
 972                                 lookup[0]->setLabel(saveLabel);
 973                                 break;
 974                         }
 975
 976                         lastLabel = lookup[0]->getLabel();
 977
 978                         //get next line to process
 979                         //prevent memory leak
 980                         for (int i = 0; i < lookup.size(); i++) {  delete lookup[i];  }
 981                         lookup = input->getSharedRAbundVectors();
 982                 }
 983
 984
 985                 if (m->control_pressed) { delete input; return 0;  }
 986
 987                 //output error messages about any remaining user labels
 988                 set<string>::iterator it;
 989                 bool needToRun = false;
 990                 for (it = userLabels.begin(); it != userLabels.end(); it++) {
 991                         m->mothurOut("Your file does not include the label " + *it);
 992                         if (processedLabels.count(lastLabel) != 1) {
 993                                 m->mothurOut(". I will use " + lastLabel + "."); m->mothurOutEndLine();
 994                                 needToRun = true;
 995                         }else {
 996                                 m->mothurOut(". Please refer to " + lastLabel + "."); m->mothurOutEndLine();
 997                         }
 998                 }
 999
1000                 //run last label if you need to
1001                 if (needToRun == true)  {
1002                         for (int i = 0; i < lookup.size(); i++) {  if (lookup[i] != NULL) {     delete lookup[i];       } }
1003                         lookup = input->getSharedRAbundVectors(lastLabel);
1004                 }
1005
1006                 delete input;
1007                 return 0;
1008         }
1009         catch(exception& e) {
1010                 m->errorOut(e, "IndicatorCommand", "getShared");
1011                 exit(1);
1012         }
1013 }
1014 //**********************************************************************************************************************
1015 int IndicatorCommand::getSharedFloat(){
1016         try {
1017                 InputData* input = new InputData(relabundfile, "relabund");
1018                 lookupFloat = input->getSharedRAbundFloatVectors();
1019                 string lastLabel = lookupFloat[0]->getLabel();
1020
1021                 if (label == "") { label = lastLabel; delete input; return 0; }
1022
1023                 //if the users enters label "0.06" and there is no "0.06" in their file use the next lowest label.
1024                 set<string> labels; labels.insert(label);
1025                 set<string> processedLabels;
1026                 set<string> userLabels = labels;
1027
1028                 //as long as you are not at the end of the file or done wih the lines you want
1029                 while((lookupFloat[0] != NULL) && (userLabels.size() != 0)) {
1030
1031                         if (m->control_pressed) {  delete input; return 0;  }
1032
1033                         if(labels.count(lookupFloat[0]->getLabel()) == 1){
1034                                 processedLabels.insert(lookupFloat[0]->getLabel());
1035                                 userLabels.erase(lookupFloat[0]->getLabel());
1036                                 break;
1037                         }
1038
1039                         if ((m->anyLabelsToProcess(lookupFloat[0]->getLabel(), userLabels, "") == true) && (processedLabels.count(lastLabel) != 1)) {
1040                                 string saveLabel = lookupFloat[0]->getLabel();
1041
1042                                 for (int i = 0; i < lookupFloat.size(); i++) {  delete lookupFloat[i];  }
1043                                 lookupFloat = input->getSharedRAbundFloatVectors(lastLabel);
1044
1045                                 processedLabels.insert(lookupFloat[0]->getLabel());
1046                                 userLabels.erase(lookupFloat[0]->getLabel());
1047
1048                                 //restore real lastlabel to save below
1049                                 lookupFloat[0]->setLabel(saveLabel);
1050                                 break;
1051                         }
1052
1053                         lastLabel = lookupFloat[0]->getLabel();
1054
1055                         //get next line to process
1056                         //prevent memory leak
1057                         for (int i = 0; i < lookupFloat.size(); i++) {  delete lookupFloat[i];  }
1058                         lookupFloat = input->getSharedRAbundFloatVectors();
1059                 }
1060
1061
1062                 if (m->control_pressed) { delete input; return 0;  }
1063
1064                 //output error messages about any remaining user labels
1065                 set<string>::iterator it;
1066                 bool needToRun = false;
1067                 for (it = userLabels.begin(); it != userLabels.end(); it++) {
1068                         m->mothurOut("Your file does not include the label " + *it);
1069                         if (processedLabels.count(lastLabel) != 1) {
1070                                 m->mothurOut(". I will use " + lastLabel + "."); m->mothurOutEndLine();
1071                                 needToRun = true;
1072                         }else {
1073                                 m->mothurOut(". Please refer to " + lastLabel + "."); m->mothurOutEndLine();
1074                         }
1075                 }
1076
1077                 //run last label if you need to
1078                 if (needToRun == true)  {
1079                         for (int i = 0; i < lookupFloat.size(); i++) {  if (lookupFloat[i] != NULL) {   delete lookupFloat[i];  } }
1080                         lookupFloat = input->getSharedRAbundFloatVectors(lastLabel);
1081                 }
1082
1083                 delete input;
1084                 return 0;
1085         }
1086         catch(exception& e) {
1087                 m->errorOut(e, "IndicatorCommand", "getShared");
1088         exit(1);
1089         }
1090 }
1091 //**********************************************************************************************************************
1092 vector<float> IndicatorCommand::driver(vector< vector<SharedRAbundFloatVector*> >& groupings, map< vector<int>, vector<int> > groupingsMap, int num, vector<float> indicatorValues, int numIters){
1093         try {
1094                 vector<float> pvalues;
1095                 pvalues.resize(indicatorValues.size(), 0);
1096
1097                 for(int i=0;i<numIters;i++){
1098                         if (m->control_pressed) { break; }
1099                         groupingsMap = randomizeGroupings(groupings, num);
1100                         vector<float> randomIndicatorValues = getValues(groupings, groupingsMap);
1101
1102                         for (int j = 0; j < indicatorValues.size(); j++) {
1103                                 if (randomIndicatorValues[j] >= indicatorValues[j]) { pvalues[j]++; }
1104                         }
1105                 }
1106
1107                 return pvalues;
1108
1109         }catch(exception& e) {
1110                 m->errorOut(e, "IndicatorCommand", "driver");
1111                 exit(1);
1112         }
1113 }
1114 //**********************************************************************************************************************
1115 vector<float> IndicatorCommand::getPValues(vector< vector<SharedRAbundFloatVector*> >& groupings, map< vector<int>, vector<int> > groupingsMap, int num, vector<float> indicatorValues){
1116         try {
1117                 vector<float> pvalues;
1118
1119 #if defined (__APPLE__) || (__MACH__) || (linux) || (__linux)
1120                 if(processors == 1){
1121                         pvalues = driver(groupings, groupingsMap, num, indicatorValues, iters);
1122                         for (int i = 0; i < pvalues.size(); i++) { pvalues[i] /= (double)iters; }
1123                 }else{
1124
1125                         //divide iters between processors
1126                         int numItersPerProcessor = iters / processors;
1127
1128                         vector<int> processIDS;
1129                         int process = 1;
1130                         int num = 0;
1131
1132                         //loop through and create all the processes you want
1133                         while (process != processors) {
1134                                 int pid = fork();
1135
1136                                 if (pid > 0) {
1137                                         processIDS.push_back(pid);  //create map from line number to pid so you can append files in correct order later
1138                                         process++;
1139                                 }else if (pid == 0){
1140                                         pvalues = driver(groupings, groupingsMap, num, indicatorValues, numItersPerProcessor);
1141
1142                                         //pass pvalues to parent
1143                                         ofstream out;
1144                                         string tempFile = toString(getpid()) + ".pvalues.temp";
1145                                         m->openOutputFile(tempFile, out);
1146
1147                                         //pass values
1148                                         for (int i = 0; i < pvalues.size(); i++) {
1149                                                 out << pvalues[i] << '\t';
1150                                         }
1151                                         out << endl;
1152
1153                                         out.close();
1154
1155                                         exit(0);
1156                                 }else {
1157                                         m->mothurOut("[ERROR]: unable to spawn the necessary processes."); m->mothurOutEndLine();
1158                                         for (int i = 0; i < processIDS.size(); i++) { kill (processIDS[i], SIGINT); }
1159                                         exit(0);
1160                                 }
1161                         }
1162
1163                         //do my part
1164                         //special case for last processor in case it doesn't divide evenly
1165                         numItersPerProcessor = iters - ((processors-1) * numItersPerProcessor);
1166                         pvalues = driver(groupings, groupingsMap, num, indicatorValues, numItersPerProcessor);
1167
1168                         //force parent to wait until all the processes are done
1169                         for (int i=0;i<processIDS.size();i++) {
1170                                 int temp = processIDS[i];
1171                                 wait(&temp);
1172                         }
1173
1174                         //combine results
1175                         for (int i = 0; i < processIDS.size(); i++) {
1176                                 ifstream in;
1177                                 string tempFile =  toString(processIDS[i]) + ".pvalues.temp";
1178                                 m->openInputFile(tempFile, in);
1179
1180                                 ////// to do ///////////
1181                                 int numTemp; numTemp = 0;
1182                                 for (int j = 0; j < pvalues.size(); j++) {
1183                                         in >> numTemp; m->gobble(in);
1184                                         pvalues[j] += numTemp;
1185                                 }
1186                                 in.close(); remove(tempFile.c_str());
1187                         }
1188                         for (int i = 0; i < pvalues.size(); i++) { pvalues[i] /= (double)iters; }
1189                 }
1190 #else
1191                 pvalues = driver(groupings, groupingsMap, num, indicatorValues, iters);
1192                 for (int i = 0; i < pvalues.size(); i++) { pvalues[i] /= (double)iters; }
1193 #endif
1194
1195                 return pvalues;
1196         }
1197         catch(exception& e) {
1198                 m->errorOut(e, "IndicatorCommand", "getPValues");
1199                 exit(1);
1200         }
1201 }
1202
1203 //**********************************************************************************************************************
1204 //same as above, just data type difference
1205 vector<float> IndicatorCommand::driver(vector< vector<SharedRAbundVector*> >& groupings, map< vector<int>, vector<int> > groupingsMap, int num, vector<float> indicatorValues, int numIters){
1206         try {
1207                 vector<float> pvalues;
1208                 pvalues.resize(indicatorValues.size(), 0);
1209
1210                 for(int i=0;i<numIters;i++){
1211                         if (m->control_pressed) { break; }
1212                         groupingsMap = randomizeGroupings(groupings, num);
1213                         vector<float> randomIndicatorValues = getValues(groupings, groupingsMap);
1214
1215                         for (int j = 0; j < indicatorValues.size(); j++) {
1216                                 if (randomIndicatorValues[j] >= indicatorValues[j]) { pvalues[j]++; }
1217                         }
1218                 }
1219
1220                 return pvalues;
1221
1222         }catch(exception& e) {
1223                 m->errorOut(e, "IndicatorCommand", "driver");
1224                 exit(1);
1225         }
1226 }
1227 //**********************************************************************************************************************
1228 //same as above, just data type difference
1229 vector<float> IndicatorCommand::getPValues(vector< vector<SharedRAbundVector*> >& groupings, map< vector<int>, vector<int> > groupingsMap, int num, vector<float> indicatorValues){
1230         try {
1231                 vector<float> pvalues;
1232
1233 #if defined (__APPLE__) || (__MACH__) || (linux) || (__linux)
1234                 if(processors == 1){
1235                         pvalues = driver(groupings, groupingsMap, num, indicatorValues, iters);
1236                         for (int i = 0; i < pvalues.size(); i++) { pvalues[i] /= (double)iters; }
1237                 }else{
1238
1239                         //divide iters between processors
1240                         int numItersPerProcessor = iters / processors;
1241
1242                         vector<int> processIDS;
1243                         int process = 1;
1244
1245                         //loop through and create all the processes you want
1246                         while (process != processors) {
1247                                 int pid = fork();
1248
1249                                 if (pid > 0) {
1250                                         processIDS.push_back(pid);  //create map from line number to pid so you can append files in correct order later
1251                                         process++;
1252                                 }else if (pid == 0){
1253                                         pvalues = driver(groupings, groupingsMap, num, indicatorValues, numItersPerProcessor);
1254
1255                                         //pass pvalues to parent
1256                                         ofstream out;
1257                                         string tempFile = toString(getpid()) + ".pvalues.temp";
1258                                         m->openOutputFile(tempFile, out);
1259
1260                                         //pass values
1261                                         for (int i = 0; i < pvalues.size(); i++) {
1262                                                 out << pvalues[i] << '\t';
1263                                         }
1264                                         out << endl;
1265
1266                                         out.close();
1267
1268                                         exit(0);
1269                                 }else {
1270                                         m->mothurOut("[ERROR]: unable to spawn the necessary processes."); m->mothurOutEndLine();
1271                                         for (int i = 0; i < processIDS.size(); i++) { kill (processIDS[i], SIGINT); }
1272                                         exit(0);
1273                                 }
1274                         }
1275
1276                         //do my part
1277                         //special case for last processor in case it doesn't divide evenly
1278                         numItersPerProcessor = iters - ((processors-1) * numItersPerProcessor);
1279                         pvalues = driver(groupings, groupingsMap, num, indicatorValues, numItersPerProcessor);
1280
1281                         //force parent to wait until all the processes are done
1282                         for (int i=0;i<processIDS.size();i++) {
1283                                 int temp = processIDS[i];
1284                                 wait(&temp);
1285                         }
1286
1287                         //combine results
1288                         for (int i = 0; i < processIDS.size(); i++) {
1289                                 ifstream in;
1290                                 string tempFile =  toString(processIDS[i]) + ".pvalues.temp";
1291                                 m->openInputFile(tempFile, in);
1292
1293                                 ////// to do ///////////
1294                                 int numTemp; numTemp = 0;
1295                                 for (int j = 0; j < pvalues.size(); j++) {
1296                                         in >> numTemp; m->gobble(in);
1297                                         pvalues[j] += numTemp;
1298                                 }
1299                                 in.close(); remove(tempFile.c_str());
1300                         }
1301                         for (int i = 0; i < pvalues.size(); i++) { pvalues[i] /= (double)iters; }
1302                 }
1303 #else
1304                 pvalues = driver(groupings, groupingsMap, num, indicatorValues, iters);
1305                 for (int i = 0; i < pvalues.size(); i++) { pvalues[i] /= (double)iters; }
1306 #endif
1307
1308                 return pvalues;
1309         }
1310         catch(exception& e) {
1311                 m->errorOut(e, "IndicatorCommand", "getPValues");
1312                 exit(1);
1313         }
1314 }
1315 //**********************************************************************************************************************
1316 //swap groups between groupings, in essence randomizing the second column of the design file
1317 map< vector<int>, vector<int> > IndicatorCommand::randomizeGroupings(vector< vector<SharedRAbundVector*> >& groupings, int numLookupGroups){
1318         try {
1319
1320                 map< vector<int>, vector<int> > randomGroupings;
1321
1322                 for (int i = 0; i < numLookupGroups; i++) {
1323                         if (m->control_pressed) {break;}
1324
1325                         //get random groups to swap to switch with
1326                         //generate random int between 0 and groupings.size()-1
1327                         int z = m->getRandomIndex(groupings.size()-1);
1328                         int x = m->getRandomIndex(groupings.size()-1);
1329                         int a = m->getRandomIndex(groupings[z].size()-1);
1330                         int b = m->getRandomIndex(groupings[x].size()-1);
1331                         //cout << i << '\t' << z << '\t' << x << '\t' << a << '\t' << b << endl;
1332                         //if ((z < 0) || (z > 1) || x<0 || x>1 || a <0 || a>groupings[z].size()-1 || b<0 || b>groupings[x].size()-1) { cout << "probelm" << i << '\t' << z << '\t' << x << '\t' << a << '\t' << b << endl;      }
1333
1334                         vector<int> from;
1335                         vector<int> to;
1336
1337                         from.push_back(z); from.push_back(a);
1338                         to.push_back(x); to.push_back(b);
1339
1340                         randomGroupings[from] = to;
1341                 }
1342                 //cout << "done" << endl;
1343                 return randomGroupings;
1344         }
1345         catch(exception& e) {
1346                 m->errorOut(e, "IndicatorCommand", "randomizeGroupings");
1347                 exit(1);
1348         }
1349 }
1350 //**********************************************************************************************************************
1351 //swap groups between groupings, in essence randomizing the second column of the design file
1352 map< vector<int>, vector<int> > IndicatorCommand::randomizeGroupings(vector< vector<SharedRAbundFloatVector*> >& groupings, int numLookupGroups){
1353         try {
1354
1355                 map< vector<int>, vector<int> > randomGroupings;
1356
1357                 for (int i = 0; i < numLookupGroups; i++) {
1358
1359                         //get random groups to swap to switch with
1360                         //generate random int between 0 and groupings.size()-1
1361                         int z = m->getRandomIndex(groupings.size()-1);
1362                         int x = m->getRandomIndex(groupings.size()-1);
1363                         int a = m->getRandomIndex(groupings[z].size()-1);
1364                         int b = m->getRandomIndex(groupings[x].size()-1);
1365                         //cout << i << '\t' << z << '\t' << x << '\t' << a << '\t' << b << endl;
1366
1367                         vector<int> from;
1368                         vector<int> to;
1369
1370                         from.push_back(z); from.push_back(a);
1371                         to.push_back(x); to.push_back(b);
1372
1373                         randomGroupings[from] = to;
1374                 }
1375
1376                 return randomGroupings;
1377         }
1378         catch(exception& e) {
1379                 m->errorOut(e, "IndicatorCommand", "randomizeGroupings");
1380                 exit(1);
1381         }
1382 }
1383
1384 /*****************************************************************/
1385
1386