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