5 * Created by westcott on 11/12/10.
6 * Copyright 2010 Schloss Lab. All rights reserved.
10 #include "indicatorcommand.h"
11 #include "sharedutilities.h"
14 //**********************************************************************************************************************
15 vector<string> IndicatorCommand::setParameters(){
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);
28 vector<string> myArray;
29 for (int i = 0; i < parameters.size(); i++) { myArray.push_back(parameters[i].name); }
33 m->errorOut(e, "IndicatorCommand", "setParameters");
37 //**********************************************************************************************************************
38 string IndicatorCommand::getHelpString(){
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";
56 m->errorOut(e, "IndicatorCommand", "getHelpString");
61 //**********************************************************************************************************************
62 IndicatorCommand::IndicatorCommand(){
64 abort = true; calledHelp = true;
66 vector<string> tempOutNames;
67 outputTypes["tree"] = tempOutNames;
68 outputTypes["summary"] = tempOutNames;
71 m->errorOut(e, "IndicatorCommand", "IndicatorCommand");
75 //**********************************************************************************************************************
76 IndicatorCommand::IndicatorCommand(string option) {
78 abort = false; calledHelp = false;
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;}
85 vector<string> myArray = setParameters();
87 OptionParser parser(option);
88 map<string, string> parameters = parser.getParameters();
90 ValidParameters validParameter;
91 map<string, string>::iterator it;
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; }
101 m->Treenames.clear();
104 vector<string> tempOutNames;
105 outputTypes["tree"] = tempOutNames;
106 outputTypes["summary"] = tempOutNames;
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 = ""; }
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; }
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; }
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; }
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; }
146 outputDir = validParameter.validFile(parameters, "outputdir", false); if (outputDir == "not found"){ outputDir = ""; }
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); }
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); }
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); }
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); }
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);
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=""; }
177 string temp = validParameter.validFile(parameters, "iters", false); if (temp == "not found") { temp = "1000"; }
178 m->mothurConvert(temp, iters);
180 temp = validParameter.validFile(parameters, "processors", false); if (temp == "not found"){ temp = m->getProcessors(); }
181 m->setProcessors(temp);
182 m->mothurConvert(temp, processors);
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(); }
190 relabundfile = m->getRelAbundFile();
191 if (relabundfile != "") { inputFileName = relabundfile; m->mothurOut("Using " + relabundfile + " as input file for the relabund parameter."); m->mothurOutEndLine(); }
193 m->mothurOut("No valid current files. You must provide a shared or relabund."); m->mothurOutEndLine();
200 if ((designfile == "") && (treefile == "")) {
201 treefile = m->getTreeFile();
202 if (treefile != "") { m->mothurOut("Using " + treefile + " as input file for the tree parameter."); m->mothurOutEndLine(); }
204 designfile = m->getDesignFile();
205 if (designfile != "") { m->mothurOut("Using " + designfile + " as input file for the design parameter."); m->mothurOutEndLine(); }
207 m->mothurOut("[ERROR]: You must provide either a tree or design file."); m->mothurOutEndLine(); abort = true;
212 if ((relabundfile != "") && (sharedfile != "")) { m->mothurOut("[ERROR]: You may not use both a shared and relabund file."); m->mothurOutEndLine(); abort = true; }
217 catch(exception& e) {
218 m->errorOut(e, "IndicatorCommand", "IndicatorCommand");
222 //**********************************************************************************************************************
224 int IndicatorCommand::execute(){
227 if (abort == true) { if (calledHelp) { return 0; } return 2; }
229 cout.setf(ios::fixed, ios::floatfield); cout.setf(ios::showpoint);
231 int start = time(NULL);
233 //read designfile if given and set up groups for read of sharedfiles
234 if (designfile != "") {
235 designMap = new GroupMap(designfile);
236 designMap->readDesignMap();
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);
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);
250 /***************************************************/
251 // use smart distancing to get right sharedRabund //
252 /***************************************************/
253 if (sharedfile != "") {
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; }
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; }
263 //reset groups if needed
264 if (designfile != "") { m->setGroups(Groups); }
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;
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();
283 treeMap->addSeq(m->Treenames[i], "Group1");
285 vector<string> myGroups; myGroups.push_back(m->Treenames[i]);
286 vector<string> myNames = designMap->getNamesSeqs(myGroups);
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();
294 treeMap->addSeq(m->Treenames[i], "Group1");
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; }
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]; } }
308 read = new ReadNewickTree(treefile);
309 int readOk = read->read(treeMap);
311 if (readOk != 0) { m->mothurOut("Read Terminated."); m->mothurOutEndLine(); delete treeMap; delete read; return 0; }
313 vector<Tree*> T = read->getTrees();
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;
324 T[0]->assembleTree();
326 /***************************************************/
327 // create ouptut tree - respecting pickedGroups //
328 /***************************************************/
329 Tree* outputTree = new Tree(m->getNumGroups(), treeMap);
331 outputTree->getSubTree(T[0], m->getGroups());
332 outputTree->assembleTree();
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]; }
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;
344 /***************************************************/
345 // get indicator species values //
346 /***************************************************/
347 GetIndicatorSpecies(outputTree);
348 delete outputTree; delete treeMap;
350 }else { //run with design file only
351 //get indicator species
352 GetIndicatorSpecies();
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]; } }
359 if (m->control_pressed) { for (int i = 0; i < outputNames.size(); i++) { m->mothurRemove(outputNames[i]); } return 0; }
361 //set tree file as new current treefile
362 if (treefile != "") {
364 itTypes = outputTypes.find("tree");
365 if (itTypes != outputTypes.end()) {
366 if ((itTypes->second).size() != 0) { current = (itTypes->second)[0]; m->setTreeFile(current); }
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();
378 catch(exception& e) {
379 m->errorOut(e, "IndicatorCommand", "execute");
383 //**********************************************************************************************************************
384 //divide shared or relabund file by groupings in the design file
385 //report all otu values to file
386 int IndicatorCommand::GetIndicatorSpecies(){
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);
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");
399 if (sharedfile != "") { numBins = lookup[0]->getNumBins(); }
400 else { numBins = lookupFloat[0]->getNumBins(); }
402 if (m->control_pressed) { out.close(); return 0; }
404 /*****************************************************/
405 //create vectors containing rabund info //
406 /*****************************************************/
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.
412 if (sharedfile != "") {
413 vector< vector<SharedRAbundVector*> > groupings;
414 set<string> groupsAlreadyAdded;
415 vector<SharedRAbundVector*> subset;
418 for (int i = 0; i < (designMap->getNamesOfGroups()).size(); i++) {
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());
427 if (subset.size() != 0) { groupings.push_back(subset); }
431 if (groupsAlreadyAdded.size() != lookup.size()) { m->mothurOut("[ERROR]: could not make proper groupings."); m->mothurOutEndLine(); }
433 indicatorValues = getValues(groupings, randomGroupingsMap);
435 pValues = getPValues(groupings, randomGroupingsMap, lookup.size(), indicatorValues);
437 vector< vector<SharedRAbundFloatVector*> > groupings;
438 set<string> groupsAlreadyAdded;
439 vector<SharedRAbundFloatVector*> subset;
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());
450 if (subset.size() != 0) { groupings.push_back(subset); }
454 if (groupsAlreadyAdded.size() != lookupFloat.size()) { m->mothurOut("[ERROR]: could not make proper groupings."); m->mothurOutEndLine(); }
456 indicatorValues = getValues(groupings, randomGroupingsMap);
458 pValues = getPValues(groupings, randomGroupingsMap, lookupFloat.size(), indicatorValues);
461 if (m->control_pressed) { out.close(); return 0; }
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++) {
470 if (m->control_pressed) { out.close(); return 0; }
472 out << (j+1) << '\t' << indicatorValues[j] << '\t';
474 if (pValues[j] > (1/(float)iters)) { out << pValues[j] << endl; }
475 else { out << "<" << (1/(float)iters) << endl; }
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();
491 catch(exception& e) {
492 m->errorOut(e, "IndicatorCommand", "GetIndicatorSpecies");
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){
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);
509 m->openOutputFile(outputFileName, out);
510 out.setf(ios::fixed, ios::floatfield); out.setf(ios::showpoint);
513 if (sharedfile != "") { numBins = lookup[0]->getNumBins(); }
514 else { numBins = lookupFloat[0]->getNumBins(); }
518 for (int i = 0; i < numBins; i++) { out << "OTU" << (i+1) << "_IndValue" << '\t' << "pValue" << '\t'; }
521 m->mothurOutEndLine(); m->mothurOut("Node\tSpecies\tIndicatorValue\tpValue\n");
523 string treeOutputDir = outputDir;
524 if (outputDir == "") { treeOutputDir += m->hasPath(treefile); }
525 string outputTreeFileName = treeOutputDir + m->getRootName(m->getSimpleName(treefile)) + "indicator.tre";
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; }
534 nodeToDescendants[i] = getDescendantList(T, i, nodeToDescendants, descendantNodes);
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);
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; }
546 /*****************************************************/
547 //create vectors containing rabund info //
548 /*****************************************************/
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.
554 if (sharedfile != "") {
555 vector< vector<SharedRAbundVector*> > groupings;
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.
562 set<string> groupsAlreadyAdded;
563 //create a grouping with my grouping
564 vector<SharedRAbundVector*> subset;
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());
574 if (count == doneCount) { break; } //quit once you get the rabunds for this grouping
576 if (subset.size() != 0) { groupings.push_back(subset); }
579 for (int j = (T->getNumNodes()-1); j >= 0; j--) {
582 if ((descendantNodes[i].count(j) == 0) && (distToRoot[j] >= distToRoot[i])) {
583 vector<SharedRAbundVector*> subset;
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());
593 if (count == doneCount) { break; } //quit once you get the rabunds for this grouping
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); }
601 if (groupsAlreadyAdded.size() != lookup.size()) { m->mothurOut("[ERROR]: could not make proper groupings."); m->mothurOutEndLine(); }
603 indicatorValues = getValues(groupings, randomGroupingsMap);
605 pValues = getPValues(groupings, randomGroupingsMap, lookup.size(), indicatorValues);
607 vector< vector<SharedRAbundFloatVector*> > groupings;
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.
614 set<string> groupsAlreadyAdded;
615 //create a grouping with my grouping
616 vector<SharedRAbundFloatVector*> subset;
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());
626 if (count == doneCount) { break; } //quit once you get the rabunds for this grouping
628 if (subset.size() != 0) { groupings.push_back(subset); }
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;
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());
642 if (count == doneCount) { break; } //quit once you get the rabunds for this grouping
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); }
650 if (groupsAlreadyAdded.size() != lookupFloat.size()) { m->mothurOut("[ERROR]: could not make proper groupings."); m->mothurOutEndLine(); }
652 indicatorValues = getValues(groupings, randomGroupingsMap);
654 pValues = getPValues(groupings, randomGroupingsMap, lookupFloat.size(), indicatorValues);
657 if (m->control_pressed) { out.close(); return 0; }
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++) {
666 if (m->control_pressed) { out.close(); return 0; }
668 if (pValues[j] < (1/(float)iters)) {
669 out << indicatorValues[j] << '\t' << '<' << (1/(float)iters) << '\t';
671 out << indicatorValues[j] << '\t' << pValues[j] << '\t';
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();
685 T->tree[i].setLabel((i+1));
690 m->openOutputFile(outputTreeFileName, outTree);
691 outputNames.push_back(outputTreeFileName); outputTypes["tree"].push_back(outputTreeFileName);
693 T->print(outTree, "both");
698 catch(exception& e) {
699 m->errorOut(e, "IndicatorCommand", "GetIndicatorSpecies");
703 //**********************************************************************************************************************
704 vector<float> IndicatorCommand::getValues(vector< vector<SharedRAbundFloatVector*> >& groupings, map< vector<int>, vector<int> > groupingsMap){
706 vector<float> values;
707 map< vector<int>, vector<int> >::iterator it;
710 for (int i = 0; i < groupings[0][0]->getNumBins(); i++) {
712 if (m->control_pressed) { return values; }
715 float AijDenominator = 0.0;
718 //get overall abundance of each grouping
719 for (int j = 0; j < groupings.size(); j++) {
721 float totalAbund = 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);
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++; }
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++; }
738 float Aij = (totalAbund / (float) groupings[j].size());
739 terms.push_back(Aij);
741 //percentage of sites represented
742 Bij.push_back(numNotZero / (float) groupings[j].size());
744 AijDenominator += Aij;
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;
753 if (thisValue > maxIndVal) { maxIndVal = thisValue; }
756 values.push_back(maxIndVal);
761 catch(exception& e) {
762 m->errorOut(e, "IndicatorCommand", "getValues");
766 //**********************************************************************************************************************
767 //same as above, just data type difference
768 vector<float> IndicatorCommand::getValues(vector< vector<SharedRAbundVector*> >& groupings, map< vector<int>, vector<int> > groupingsMap){
770 vector<float> values;
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;
778 map< vector<int>, vector<int> >::iterator it;
781 for (int i = 0; i < groupings[0][0]->getNumBins(); i++) {
783 float AijDenominator = 0.0;
785 //get overall abundance of each grouping
786 for (int j = 0; j < groupings.size(); j++) {
788 int totalAbund = 0.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);
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++; }
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++; }
805 float Aij = (totalAbund / (float) groupings[j].size());
806 terms.push_back(Aij);
808 //percentage of sites represented
809 Bij.push_back(numNotZero / (float) groupings[j].size());
811 AijDenominator += Aij;
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;
820 if (thisValue > maxIndVal) { maxIndVal = thisValue; }
823 values.push_back(maxIndVal);
828 catch(exception& e) {
829 m->errorOut(e, "IndicatorCommand", "getValues");
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){
838 map<int, float> dists;
840 bool hasBranchLengths = false;
841 for (int i = 0; i < T->getNumNodes(); i++) {
842 if (T->tree[i].getBranchLength() > 0.0) { hasBranchLengths = true; break; }
845 //set branchlengths if needed
846 if (!hasBranchLengths) {
847 for (int i = 0; i < T->getNumNodes(); i++) {
849 int lc = T->tree[i].getLChild();
850 int rc = T->tree[i].getRChild();
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);
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];
861 float greater = ldist;
862 if (rdist > greater) { greater = rdist; dists[i] = ldist + 1.0;}
863 else { dists[i] = rdist + 1.0; }
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));
875 for (int i = 0; i < T->getNumNodes(); i++) {
880 while(T->tree[index].getParent() != -1){
881 if (T->tree[index].getBranchLength() != -1) {
882 sum += abs(T->tree[index].getBranchLength());
884 index = T->tree[index].getParent();
892 catch(exception& e) {
893 m->errorOut(e, "IndicatorCommand", "getLengthToLeaf");
897 //**********************************************************************************************************************
898 set<string> IndicatorCommand::getDescendantList(Tree*& T, int i, map<int, set<string> > descendants, map<int, set<int> >& nodes){
902 set<string>::iterator it;
904 int lc = T->tree[i].getLChild();
905 int rc = T->tree[i].getRChild();
907 if (lc == -1) { //you are a leaf your only descendant is yourself
908 set<int> temp; temp.insert(i);
911 if (designfile == "") {
912 names.insert(T->tree[i].getName());
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]);
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++) {
927 for (set<int>::iterator itNum = nodes[rc].begin(); itNum != nodes[rc].end(); itNum++) {
928 nodes[i].insert(*itNum);
930 //you are your own descendant
936 catch(exception& e) {
937 m->errorOut(e, "IndicatorCommand", "getDescendantList");
941 //**********************************************************************************************************************
942 int IndicatorCommand::getShared(){
944 InputData* input = new InputData(sharedfile, "sharedfile");
945 lookup = input->getSharedRAbundVectors();
946 string lastLabel = lookup[0]->getLabel();
948 if (label == "") { label = lastLabel; delete input; return 0; }
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;
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; }
959 if(labels.count(lookup[0]->getLabel()) == 1){
960 processedLabels.insert(lookup[0]->getLabel());
961 userLabels.erase(lookup[0]->getLabel());
965 if ((m->anyLabelsToProcess(lookup[0]->getLabel(), userLabels, "") == true) && (processedLabels.count(lastLabel) != 1)) {
966 string saveLabel = lookup[0]->getLabel();
968 for (int i = 0; i < lookup.size(); i++) { delete lookup[i]; }
969 lookup = input->getSharedRAbundVectors(lastLabel);
971 processedLabels.insert(lookup[0]->getLabel());
972 userLabels.erase(lookup[0]->getLabel());
974 //restore real lastlabel to save below
975 lookup[0]->setLabel(saveLabel);
979 lastLabel = lookup[0]->getLabel();
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();
988 if (m->control_pressed) { delete input; return 0; }
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();
999 m->mothurOut(". Please refer to " + lastLabel + "."); m->mothurOutEndLine();
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);
1012 catch(exception& e) {
1013 m->errorOut(e, "IndicatorCommand", "getShared");
1017 //**********************************************************************************************************************
1018 int IndicatorCommand::getSharedFloat(){
1020 InputData* input = new InputData(relabundfile, "relabund");
1021 lookupFloat = input->getSharedRAbundFloatVectors();
1022 string lastLabel = lookupFloat[0]->getLabel();
1024 if (label == "") { label = lastLabel; delete input; return 0; }
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;
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)) {
1034 if (m->control_pressed) { delete input; return 0; }
1036 if(labels.count(lookupFloat[0]->getLabel()) == 1){
1037 processedLabels.insert(lookupFloat[0]->getLabel());
1038 userLabels.erase(lookupFloat[0]->getLabel());
1042 if ((m->anyLabelsToProcess(lookupFloat[0]->getLabel(), userLabels, "") == true) && (processedLabels.count(lastLabel) != 1)) {
1043 string saveLabel = lookupFloat[0]->getLabel();
1045 for (int i = 0; i < lookupFloat.size(); i++) { delete lookupFloat[i]; }
1046 lookupFloat = input->getSharedRAbundFloatVectors(lastLabel);
1048 processedLabels.insert(lookupFloat[0]->getLabel());
1049 userLabels.erase(lookupFloat[0]->getLabel());
1051 //restore real lastlabel to save below
1052 lookupFloat[0]->setLabel(saveLabel);
1056 lastLabel = lookupFloat[0]->getLabel();
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();
1065 if (m->control_pressed) { delete input; return 0; }
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();
1076 m->mothurOut(". Please refer to " + lastLabel + "."); m->mothurOutEndLine();
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);
1089 catch(exception& e) {
1090 m->errorOut(e, "IndicatorCommand", "getShared");
1094 //**********************************************************************************************************************
1095 vector<float> IndicatorCommand::driver(vector< vector<SharedRAbundFloatVector*> >& groupings, map< vector<int>, vector<int> > groupingsMap, int num, vector<float> indicatorValues, int numIters){
1097 vector<float> pvalues;
1098 pvalues.resize(indicatorValues.size(), 0);
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);
1105 for (int j = 0; j < indicatorValues.size(); j++) {
1106 if (randomIndicatorValues[j] >= indicatorValues[j]) { pvalues[j]++; }
1112 }catch(exception& e) {
1113 m->errorOut(e, "IndicatorCommand", "driver");
1117 //**********************************************************************************************************************
1118 vector<float> IndicatorCommand::getPValues(vector< vector<SharedRAbundFloatVector*> >& groupings, map< vector<int>, vector<int> > groupingsMap, int num, vector<float> indicatorValues){
1120 vector<float> pvalues;
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; }
1128 //divide iters between processors
1129 int numItersPerProcessor = iters / processors;
1131 vector<int> processIDS;
1135 //loop through and create all the processes you want
1136 while (process != processors) {
1140 processIDS.push_back(pid); //create map from line number to pid so you can append files in correct order later
1142 }else if (pid == 0){
1143 pvalues = driver(groupings, groupingsMap, num, indicatorValues, numItersPerProcessor);
1145 //pass pvalues to parent
1147 string tempFile = toString(getpid()) + ".pvalues.temp";
1148 m->openOutputFile(tempFile, out);
1151 for (int i = 0; i < pvalues.size(); i++) {
1152 out << pvalues[i] << '\t';
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); }
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);
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];
1178 for (int i = 0; i < processIDS.size(); i++) {
1180 string tempFile = toString(processIDS[i]) + ".pvalues.temp";
1181 m->openInputFile(tempFile, in);
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;
1189 in.close(); m->mothurRemove(tempFile);
1191 for (int i = 0; i < pvalues.size(); i++) { pvalues[i] /= (double)iters; }
1194 pvalues = driver(groupings, groupingsMap, num, indicatorValues, iters);
1195 for (int i = 0; i < pvalues.size(); i++) { pvalues[i] /= (double)iters; }
1200 catch(exception& e) {
1201 m->errorOut(e, "IndicatorCommand", "getPValues");
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){
1210 vector<float> pvalues;
1211 pvalues.resize(indicatorValues.size(), 0);
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);
1218 for (int j = 0; j < indicatorValues.size(); j++) {
1219 if (randomIndicatorValues[j] >= indicatorValues[j]) { pvalues[j]++; }
1225 }catch(exception& e) {
1226 m->errorOut(e, "IndicatorCommand", "driver");
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){
1234 vector<float> pvalues;
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; }
1242 //divide iters between processors
1243 int numItersPerProcessor = iters / processors;
1245 vector<int> processIDS;
1248 //loop through and create all the processes you want
1249 while (process != processors) {
1253 processIDS.push_back(pid); //create map from line number to pid so you can append files in correct order later
1255 }else if (pid == 0){
1256 pvalues = driver(groupings, groupingsMap, num, indicatorValues, numItersPerProcessor);
1258 //pass pvalues to parent
1260 string tempFile = toString(getpid()) + ".pvalues.temp";
1261 m->openOutputFile(tempFile, out);
1264 for (int i = 0; i < pvalues.size(); i++) {
1265 out << pvalues[i] << '\t';
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); }
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);
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];
1291 for (int i = 0; i < processIDS.size(); i++) {
1293 string tempFile = toString(processIDS[i]) + ".pvalues.temp";
1294 m->openInputFile(tempFile, in);
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;
1302 in.close(); m->mothurRemove(tempFile);
1304 for (int i = 0; i < pvalues.size(); i++) { pvalues[i] /= (double)iters; }
1307 pvalues = driver(groupings, groupingsMap, num, indicatorValues, iters);
1308 for (int i = 0; i < pvalues.size(); i++) { pvalues[i] /= (double)iters; }
1313 catch(exception& e) {
1314 m->errorOut(e, "IndicatorCommand", "getPValues");
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){
1323 map< vector<int>, vector<int> > randomGroupings;
1325 for (int i = 0; i < numLookupGroups; i++) {
1326 if (m->control_pressed) {break;}
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; }
1340 from.push_back(z); from.push_back(a);
1341 to.push_back(x); to.push_back(b);
1343 randomGroupings[from] = to;
1345 //cout << "done" << endl;
1346 return randomGroupings;
1348 catch(exception& e) {
1349 m->errorOut(e, "IndicatorCommand", "randomizeGroupings");
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){
1358 map< vector<int>, vector<int> > randomGroupings;
1360 for (int i = 0; i < numLookupGroups; i++) {
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;
1373 from.push_back(z); from.push_back(a);
1374 to.push_back(x); to.push_back(b);
1376 randomGroupings[from] = to;
1379 return randomGroupings;
1381 catch(exception& e) {
1382 m->errorOut(e, "IndicatorCommand", "randomizeGroupings");
1387 /*****************************************************************/