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