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]; }
338 if (m->control_pressed) {
339 if (designfile != "") { delete designMap; }
340 if (sharedfile != "") { for (int i = 0; i < lookup.size(); i++) { delete lookup[i]; } }
341 else { for (int i = 0; i < lookupFloat.size(); i++) { delete lookupFloat[i]; } }
342 delete outputTree; delete treeMap; return 0;
345 /***************************************************/
346 // get indicator species values //
347 /***************************************************/
348 GetIndicatorSpecies(outputTree);
349 delete outputTree; delete treeMap;
351 }else { //run with design file only
352 //get indicator species
353 GetIndicatorSpecies();
356 if (designfile != "") { delete designMap; }
357 if (sharedfile != "") { for (int i = 0; i < lookup.size(); i++) { delete lookup[i]; } }
358 else { for (int i = 0; i < lookupFloat.size(); i++) { delete lookupFloat[i]; } }
360 if (m->control_pressed) { for (int i = 0; i < outputNames.size(); i++) { m->mothurRemove(outputNames[i]); } return 0; }
362 //set tree file as new current treefile
363 if (treefile != "") {
365 itTypes = outputTypes.find("tree");
366 if (itTypes != outputTypes.end()) {
367 if ((itTypes->second).size() != 0) { current = (itTypes->second)[0]; m->setTreeFile(current); }
371 m->mothurOutEndLine(); m->mothurOutEndLine(); m->mothurOut("It took " + toString(time(NULL) - start) + " secs to find the indicator species."); m->mothurOutEndLine();
372 m->mothurOutEndLine();
373 m->mothurOut("Output File Names: "); m->mothurOutEndLine();
374 for (int i = 0; i < outputNames.size(); i++) { m->mothurOut(outputNames[i]); m->mothurOutEndLine(); }
375 m->mothurOutEndLine();
379 catch(exception& e) {
380 m->errorOut(e, "IndicatorCommand", "execute");
384 //**********************************************************************************************************************
385 //divide shared or relabund file by groupings in the design file
386 //report all otu values to file
387 int IndicatorCommand::GetIndicatorSpecies(){
389 string thisOutputDir = outputDir;
390 if (outputDir == "") { thisOutputDir += m->hasPath(inputFileName); }
391 string outputFileName = thisOutputDir + m->getRootName(m->getSimpleName(inputFileName)) + "indicator.summary";
392 outputNames.push_back(outputFileName); outputTypes["summary"].push_back(outputFileName);
395 m->openOutputFile(outputFileName, out);
396 out.setf(ios::fixed, ios::floatfield); out.setf(ios::showpoint);
397 m->mothurOutEndLine(); m->mothurOut("Species\tIndicatorValue\tpValue\n");
400 if (sharedfile != "") { numBins = lookup[0]->getNumBins(); }
401 else { numBins = lookupFloat[0]->getNumBins(); }
403 if (m->control_pressed) { out.close(); return 0; }
405 /*****************************************************/
406 //create vectors containing rabund info //
407 /*****************************************************/
409 vector<float> indicatorValues; //size of numBins
410 vector<float> pValues;
411 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.
413 if (sharedfile != "") {
414 vector< vector<SharedRAbundVector*> > groupings;
415 set<string> groupsAlreadyAdded;
416 vector<SharedRAbundVector*> subset;
419 for (int i = 0; i < (designMap->getNamesOfGroups()).size(); i++) {
421 for (int k = 0; k < lookup.size(); k++) {
422 //are you from this grouping?
423 if (designMap->getGroup(lookup[k]->getGroup()) == (designMap->getNamesOfGroups())[i]) {
424 subset.push_back(lookup[k]);
425 groupsAlreadyAdded.insert(lookup[k]->getGroup());
428 if (subset.size() != 0) { groupings.push_back(subset); }
432 if (groupsAlreadyAdded.size() != lookup.size()) { m->mothurOut("[ERROR]: could not make proper groupings."); m->mothurOutEndLine(); }
434 indicatorValues = getValues(groupings, randomGroupingsMap);
436 pValues = getPValues(groupings, randomGroupingsMap, lookup.size(), indicatorValues);
438 vector< vector<SharedRAbundFloatVector*> > groupings;
439 set<string> groupsAlreadyAdded;
440 vector<SharedRAbundFloatVector*> subset;
443 for (int i = 0; i < (designMap->getNamesOfGroups()).size(); i++) {
444 for (int k = 0; k < lookupFloat.size(); k++) {
445 //are you from this grouping?
446 if (designMap->getGroup(lookupFloat[k]->getGroup()) == (designMap->getNamesOfGroups())[i]) {
447 subset.push_back(lookupFloat[k]);
448 groupsAlreadyAdded.insert(lookupFloat[k]->getGroup());
451 if (subset.size() != 0) { groupings.push_back(subset); }
455 if (groupsAlreadyAdded.size() != lookupFloat.size()) { m->mothurOut("[ERROR]: could not make proper groupings."); m->mothurOutEndLine(); }
457 indicatorValues = getValues(groupings, randomGroupingsMap);
459 pValues = getPValues(groupings, randomGroupingsMap, lookupFloat.size(), indicatorValues);
462 if (m->control_pressed) { out.close(); return 0; }
465 /******************************************************/
466 //output indicator values to table form //
467 /*****************************************************/
468 out << "OTU\tIndicator_Value\tpValue" << endl;
469 for (int j = 0; j < indicatorValues.size(); j++) {
471 if (m->control_pressed) { out.close(); return 0; }
473 out << (j+1) << '\t' << indicatorValues[j] << '\t';
475 if (pValues[j] > (1/(float)iters)) { out << pValues[j] << endl; }
476 else { out << "<" << (1/(float)iters) << endl; }
478 if (pValues[j] <= 0.05) {
479 cout << "OTU" << j+1 << '\t' << indicatorValues[j] << '\t';
480 string pValueString = "<" + toString((1/(float)iters));
481 if (pValues[j] > (1/(float)iters)) { pValueString = toString(pValues[j]); cout << pValues[j];}
482 else { cout << "<" << (1/(float)iters); }
483 m->mothurOutJustToLog("OTU" + toString(j+1) + "\t" + toString(indicatorValues[j]) + "\t" + pValueString);
484 m->mothurOutEndLine();
492 catch(exception& e) {
493 m->errorOut(e, "IndicatorCommand", "GetIndicatorSpecies");
497 //**********************************************************************************************************************
498 //traverse tree finding indicator species values for each otu at each node
499 //label node with otu number that has highest indicator value
500 //report all otu values to file
501 int IndicatorCommand::GetIndicatorSpecies(Tree*& T){
504 string thisOutputDir = outputDir;
505 if (outputDir == "") { thisOutputDir += m->hasPath(inputFileName); }
506 string outputFileName = thisOutputDir + m->getRootName(m->getSimpleName(inputFileName)) + "indicator.summary";
507 outputNames.push_back(outputFileName); outputTypes["summary"].push_back(outputFileName);
510 m->openOutputFile(outputFileName, out);
511 out.setf(ios::fixed, ios::floatfield); out.setf(ios::showpoint);
514 if (sharedfile != "") { numBins = lookup[0]->getNumBins(); }
515 else { numBins = lookupFloat[0]->getNumBins(); }
519 for (int i = 0; i < numBins; i++) { out << "OTU" << (i+1) << "_IndValue" << '\t' << "pValue" << '\t'; }
522 m->mothurOutEndLine(); m->mothurOut("Node\tSpecies\tIndicatorValue\tpValue\n");
524 string treeOutputDir = outputDir;
525 if (outputDir == "") { treeOutputDir += m->hasPath(treefile); }
526 string outputTreeFileName = treeOutputDir + m->getRootName(m->getSimpleName(treefile)) + "indicator.tre";
529 //create a map from tree node index to names of descendants, save time later to know which sharedRabund you need
530 map<int, set<string> > nodeToDescendants;
531 map<int, set<int> > descendantNodes;
532 for (int i = 0; i < T->getNumNodes(); i++) {
533 if (m->control_pressed) { return 0; }
535 nodeToDescendants[i] = getDescendantList(T, i, nodeToDescendants, descendantNodes);
538 //you need the distances to leaf to decide grouping below
539 //this will also set branch lengths if the tree does not include them
540 map<int, float> distToRoot = getDistToRoot(T);
543 for (int i = T->getNumLeaves(); i < T->getNumNodes(); i++) {
544 //cout << endl << i+1 << endl;
545 if (m->control_pressed) { out.close(); return 0; }
547 /*****************************************************/
548 //create vectors containing rabund info //
549 /*****************************************************/
551 vector<float> indicatorValues; //size of numBins
552 vector<float> pValues;
553 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.
555 if (sharedfile != "") {
556 vector< vector<SharedRAbundVector*> > groupings;
558 //get nodes that will be a valid grouping
559 //you are valid if you are not one of my descendants
560 //AND your distToRoot is >= mine
561 //AND you were not added as part of a larger grouping. Largest nodes are added first.
563 set<string> groupsAlreadyAdded;
564 //create a grouping with my grouping
565 vector<SharedRAbundVector*> subset;
567 int doneCount = nodeToDescendants[i].size();
568 for (int k = 0; k < lookup.size(); k++) {
569 //is this descendant of i
570 if ((nodeToDescendants[i].count(lookup[k]->getGroup()) != 0)) {
571 subset.push_back(lookup[k]);
572 groupsAlreadyAdded.insert(lookup[k]->getGroup());
575 if (count == doneCount) { break; } //quit once you get the rabunds for this grouping
577 if (subset.size() != 0) { groupings.push_back(subset); }
580 for (int j = (T->getNumNodes()-1); j >= 0; j--) {
583 if ((descendantNodes[i].count(j) == 0) && (distToRoot[j] >= distToRoot[i])) {
584 vector<SharedRAbundVector*> subset;
586 int doneCount = nodeToDescendants[j].size();
587 for (int k = 0; k < lookup.size(); k++) {
588 //is this descendant of j, and we didn't already add this as part of a larger grouping
589 if ((nodeToDescendants[j].count(lookup[k]->getGroup()) != 0) && (groupsAlreadyAdded.count(lookup[k]->getGroup()) == 0)) {
590 subset.push_back(lookup[k]);
591 groupsAlreadyAdded.insert(lookup[k]->getGroup());
594 if (count == doneCount) { break; } //quit once you get the rabunds for this grouping
597 //if subset.size == 0 then the node was added as part of a larger grouping
598 if (subset.size() != 0) { groupings.push_back(subset); }
602 if (groupsAlreadyAdded.size() != lookup.size()) { m->mothurOut("[ERROR]: could not make proper groupings."); m->mothurOutEndLine(); }
604 indicatorValues = getValues(groupings, randomGroupingsMap);
606 pValues = getPValues(groupings, randomGroupingsMap, lookup.size(), indicatorValues);
608 vector< vector<SharedRAbundFloatVector*> > groupings;
610 //get nodes that will be a valid grouping
611 //you are valid if you are not one of my descendants
612 //AND your distToRoot is >= mine
613 //AND you were not added as part of a larger grouping. Largest nodes are added first.
615 set<string> groupsAlreadyAdded;
616 //create a grouping with my grouping
617 vector<SharedRAbundFloatVector*> subset;
619 int doneCount = nodeToDescendants[i].size();
620 for (int k = 0; k < lookupFloat.size(); k++) {
621 //is this descendant of i
622 if ((nodeToDescendants[i].count(lookupFloat[k]->getGroup()) != 0)) {
623 subset.push_back(lookupFloat[k]);
624 groupsAlreadyAdded.insert(lookupFloat[k]->getGroup());
627 if (count == doneCount) { break; } //quit once you get the rabunds for this grouping
629 if (subset.size() != 0) { groupings.push_back(subset); }
631 for (int j = (T->getNumNodes()-1); j >= 0; j--) {
632 if ((descendantNodes[i].count(j) == 0) && (distToRoot[j] >= distToRoot[i])) {
633 vector<SharedRAbundFloatVector*> subset;
635 int doneCount = nodeToDescendants[j].size();
636 for (int k = 0; k < lookupFloat.size(); k++) {
637 //is this descendant of j, and we didn't already add this as part of a larger grouping
638 if ((nodeToDescendants[j].count(lookupFloat[k]->getGroup()) != 0) && (groupsAlreadyAdded.count(lookupFloat[k]->getGroup()) == 0)) {
639 subset.push_back(lookupFloat[k]);
640 groupsAlreadyAdded.insert(lookupFloat[k]->getGroup());
643 if (count == doneCount) { break; } //quit once you get the rabunds for this grouping
646 //if subset.size == 0 then the node was added as part of a larger grouping
647 if (subset.size() != 0) { groupings.push_back(subset); }
651 if (groupsAlreadyAdded.size() != lookupFloat.size()) { m->mothurOut("[ERROR]: could not make proper groupings."); m->mothurOutEndLine(); }
653 indicatorValues = getValues(groupings, randomGroupingsMap);
655 pValues = getPValues(groupings, randomGroupingsMap, lookupFloat.size(), indicatorValues);
658 if (m->control_pressed) { out.close(); return 0; }
661 /******************************************************/
662 //output indicator values to table form + label tree //
663 /*****************************************************/
664 out << (i+1) << '\t';
665 for (int j = 0; j < indicatorValues.size(); j++) {
667 if (m->control_pressed) { out.close(); return 0; }
669 if (pValues[j] < (1/(float)iters)) {
670 out << indicatorValues[j] << '\t' << '<' << (1/(float)iters) << '\t';
672 out << indicatorValues[j] << '\t' << pValues[j] << '\t';
675 if (pValues[j] <= 0.05) {
676 cout << i+1 << "\tOTU" << j+1 << '\t' << indicatorValues[j] << '\t';
677 string pValueString = "<" + toString((1/(float)iters));
678 if (pValues[j] > (1/(float)iters)) { pValueString = toString(pValues[j]); cout << pValues[j];}
679 else { cout << "<" << (1/(float)iters); }
680 m->mothurOutJustToLog(toString(i) + "\tOTU" + toString(j+1) + "\t" + toString(indicatorValues[j]) + "\t" + pValueString);
681 m->mothurOutEndLine();
686 T->tree[i].setLabel((i+1));
691 m->openOutputFile(outputTreeFileName, outTree);
692 outputNames.push_back(outputTreeFileName); outputTypes["tree"].push_back(outputTreeFileName);
694 T->print(outTree, "both");
699 catch(exception& e) {
700 m->errorOut(e, "IndicatorCommand", "GetIndicatorSpecies");
704 //**********************************************************************************************************************
705 vector<float> IndicatorCommand::getValues(vector< vector<SharedRAbundFloatVector*> >& groupings, map< vector<int>, vector<int> > groupingsMap){
707 vector<float> values;
708 map< vector<int>, vector<int> >::iterator it;
711 for (int i = 0; i < groupings[0][0]->getNumBins(); i++) {
713 if (m->control_pressed) { return values; }
716 float AijDenominator = 0.0;
719 //get overall abundance of each grouping
720 for (int j = 0; j < groupings.size(); j++) {
722 float totalAbund = 0;
724 for (int k = 0; k < groupings[j].size(); k++) {
725 vector<int> temp; temp.push_back(j); temp.push_back(k);
726 it = groupingsMap.find(temp);
728 if (it == groupingsMap.end()) { //this one didnt get moved
729 totalAbund += groupings[j][k]->getAbundance(i);
730 if (groupings[j][k]->getAbundance(i) != 0.0) { numNotZero++; }
732 totalAbund += groupings[(it->second)[0]][(it->second)[1]]->getAbundance(i);
733 if (groupings[(it->second)[0]][(it->second)[1]]->getAbundance(i) != 0.0) { numNotZero++; }
739 float Aij = (totalAbund / (float) groupings[j].size());
740 terms.push_back(Aij);
742 //percentage of sites represented
743 Bij.push_back(numNotZero / (float) groupings[j].size());
745 AijDenominator += Aij;
748 float maxIndVal = 0.0;
749 for (int j = 0; j < terms.size(); j++) {
750 float thisAij = (terms[j] / AijDenominator); //relative abundance
751 float thisValue = thisAij * Bij[j] * 100.0;
754 if (thisValue > maxIndVal) { maxIndVal = thisValue; }
757 values.push_back(maxIndVal);
762 catch(exception& e) {
763 m->errorOut(e, "IndicatorCommand", "getValues");
767 //**********************************************************************************************************************
768 //same as above, just data type difference
769 vector<float> IndicatorCommand::getValues(vector< vector<SharedRAbundVector*> >& groupings, map< vector<int>, vector<int> > groupingsMap){
771 vector<float> values;
773 /*for (int j = 0; j < groupings.size(); j++) {
774 cout << "grouping " << j << endl;
775 for (int k = 0; k < groupings[j].size(); k++) {
776 cout << groupings[j][k]->getGroup() << endl;
779 map< vector<int>, vector<int> >::iterator it;
782 for (int i = 0; i < groupings[0][0]->getNumBins(); i++) {
784 float AijDenominator = 0.0;
786 //get overall abundance of each grouping
787 for (int j = 0; j < groupings.size(); j++) {
789 int totalAbund = 0.0;
791 for (int k = 0; k < groupings[j].size(); k++) {
792 vector<int> temp; temp.push_back(j); temp.push_back(k);
793 it = groupingsMap.find(temp);
795 if (it == groupingsMap.end()) { //this one didnt get moved
796 totalAbund += groupings[j][k]->getAbundance(i);
797 if (groupings[j][k]->getAbundance(i) != 0.0) { numNotZero++; }
799 totalAbund += groupings[(it->second)[0]][(it->second)[1]]->getAbundance(i);
800 if (groupings[(it->second)[0]][(it->second)[1]]->getAbundance(i) != 0.0) { numNotZero++; }
806 float Aij = (totalAbund / (float) groupings[j].size());
807 terms.push_back(Aij);
809 //percentage of sites represented
810 Bij.push_back(numNotZero / (float) groupings[j].size());
812 AijDenominator += Aij;
815 float maxIndVal = 0.0;
816 for (int j = 0; j < terms.size(); j++) {
817 float thisAij = (terms[j] / AijDenominator); //relative abundance
818 float thisValue = thisAij * Bij[j] * 100.0;
821 if (thisValue > maxIndVal) { maxIndVal = thisValue; }
824 values.push_back(maxIndVal);
829 catch(exception& e) {
830 m->errorOut(e, "IndicatorCommand", "getValues");
834 //**********************************************************************************************************************
835 //you need the distances to root to decide groupings
836 //this will also set branch lengths if the tree does not include them
837 map<int, float> IndicatorCommand::getDistToRoot(Tree*& T){
839 map<int, float> dists;
841 bool hasBranchLengths = false;
842 for (int i = 0; i < T->getNumNodes(); i++) {
843 if (T->tree[i].getBranchLength() > 0.0) { hasBranchLengths = true; break; }
846 //set branchlengths if needed
847 if (!hasBranchLengths) {
848 for (int i = 0; i < T->getNumNodes(); i++) {
850 int lc = T->tree[i].getLChild();
851 int rc = T->tree[i].getRChild();
853 if (lc == -1) { // you are a leaf
854 //if you are a leaf set you priliminary length to 1.0, this may adjust later
855 T->tree[i].setBranchLength(1.0);
857 }else{ // you are an internal node
858 //look at your children's length to leaf
859 float ldist = dists[lc];
860 float rdist = dists[rc];
862 float greater = ldist;
863 if (rdist > greater) { greater = rdist; dists[i] = ldist + 1.0;}
864 else { dists[i] = rdist + 1.0; }
867 //branch length = difference + 1
868 T->tree[lc].setBranchLength((abs(ldist-greater) + 1.0));
869 T->tree[rc].setBranchLength((abs(rdist-greater) + 1.0));
876 for (int i = 0; i < T->getNumNodes(); i++) {
881 while(T->tree[index].getParent() != -1){
882 if (T->tree[index].getBranchLength() != -1) {
883 sum += abs(T->tree[index].getBranchLength());
885 index = T->tree[index].getParent();
893 catch(exception& e) {
894 m->errorOut(e, "IndicatorCommand", "getLengthToLeaf");
898 //**********************************************************************************************************************
899 set<string> IndicatorCommand::getDescendantList(Tree*& T, int i, map<int, set<string> > descendants, map<int, set<int> >& nodes){
903 set<string>::iterator it;
905 int lc = T->tree[i].getLChild();
906 int rc = T->tree[i].getRChild();
908 if (lc == -1) { //you are a leaf your only descendant is yourself
909 set<int> temp; temp.insert(i);
912 if (designfile == "") {
913 names.insert(T->tree[i].getName());
915 vector<string> myGroup; myGroup.push_back(T->tree[i].getName());
916 vector<string> myReps = designMap->getNamesSeqs(myGroup);
917 for (int k = 0; k < myReps.size(); k++) {
918 names.insert(myReps[k]);
922 }else{ //your descedants are the combination of your childrens descendants
923 names = descendants[lc];
924 nodes[i] = nodes[lc];
925 for (it = descendants[rc].begin(); it != descendants[rc].end(); it++) {
928 for (set<int>::iterator itNum = nodes[rc].begin(); itNum != nodes[rc].end(); itNum++) {
929 nodes[i].insert(*itNum);
931 //you are your own descendant
937 catch(exception& e) {
938 m->errorOut(e, "IndicatorCommand", "getDescendantList");
942 //**********************************************************************************************************************
943 int IndicatorCommand::getShared(){
945 InputData* input = new InputData(sharedfile, "sharedfile");
946 lookup = input->getSharedRAbundVectors();
947 string lastLabel = lookup[0]->getLabel();
949 if (label == "") { label = lastLabel; delete input; return 0; }
951 //if the users enters label "0.06" and there is no "0.06" in their file use the next lowest label.
952 set<string> labels; labels.insert(label);
953 set<string> processedLabels;
954 set<string> userLabels = labels;
956 //as long as you are not at the end of the file or done wih the lines you want
957 while((lookup[0] != NULL) && (userLabels.size() != 0)) {
958 if (m->control_pressed) { delete input; return 0; }
960 if(labels.count(lookup[0]->getLabel()) == 1){
961 processedLabels.insert(lookup[0]->getLabel());
962 userLabels.erase(lookup[0]->getLabel());
966 if ((m->anyLabelsToProcess(lookup[0]->getLabel(), userLabels, "") == true) && (processedLabels.count(lastLabel) != 1)) {
967 string saveLabel = lookup[0]->getLabel();
969 for (int i = 0; i < lookup.size(); i++) { delete lookup[i]; }
970 lookup = input->getSharedRAbundVectors(lastLabel);
972 processedLabels.insert(lookup[0]->getLabel());
973 userLabels.erase(lookup[0]->getLabel());
975 //restore real lastlabel to save below
976 lookup[0]->setLabel(saveLabel);
980 lastLabel = lookup[0]->getLabel();
982 //get next line to process
983 //prevent memory leak
984 for (int i = 0; i < lookup.size(); i++) { delete lookup[i]; }
985 lookup = input->getSharedRAbundVectors();
989 if (m->control_pressed) { delete input; return 0; }
991 //output error messages about any remaining user labels
992 set<string>::iterator it;
993 bool needToRun = false;
994 for (it = userLabels.begin(); it != userLabels.end(); it++) {
995 m->mothurOut("Your file does not include the label " + *it);
996 if (processedLabels.count(lastLabel) != 1) {
997 m->mothurOut(". I will use " + lastLabel + "."); m->mothurOutEndLine();
1000 m->mothurOut(". Please refer to " + lastLabel + "."); m->mothurOutEndLine();
1004 //run last label if you need to
1005 if (needToRun == true) {
1006 for (int i = 0; i < lookup.size(); i++) { if (lookup[i] != NULL) { delete lookup[i]; } }
1007 lookup = input->getSharedRAbundVectors(lastLabel);
1013 catch(exception& e) {
1014 m->errorOut(e, "IndicatorCommand", "getShared");
1018 //**********************************************************************************************************************
1019 int IndicatorCommand::getSharedFloat(){
1021 InputData* input = new InputData(relabundfile, "relabund");
1022 lookupFloat = input->getSharedRAbundFloatVectors();
1023 string lastLabel = lookupFloat[0]->getLabel();
1025 if (label == "") { label = lastLabel; delete input; return 0; }
1027 //if the users enters label "0.06" and there is no "0.06" in their file use the next lowest label.
1028 set<string> labels; labels.insert(label);
1029 set<string> processedLabels;
1030 set<string> userLabels = labels;
1032 //as long as you are not at the end of the file or done wih the lines you want
1033 while((lookupFloat[0] != NULL) && (userLabels.size() != 0)) {
1035 if (m->control_pressed) { delete input; return 0; }
1037 if(labels.count(lookupFloat[0]->getLabel()) == 1){
1038 processedLabels.insert(lookupFloat[0]->getLabel());
1039 userLabels.erase(lookupFloat[0]->getLabel());
1043 if ((m->anyLabelsToProcess(lookupFloat[0]->getLabel(), userLabels, "") == true) && (processedLabels.count(lastLabel) != 1)) {
1044 string saveLabel = lookupFloat[0]->getLabel();
1046 for (int i = 0; i < lookupFloat.size(); i++) { delete lookupFloat[i]; }
1047 lookupFloat = input->getSharedRAbundFloatVectors(lastLabel);
1049 processedLabels.insert(lookupFloat[0]->getLabel());
1050 userLabels.erase(lookupFloat[0]->getLabel());
1052 //restore real lastlabel to save below
1053 lookupFloat[0]->setLabel(saveLabel);
1057 lastLabel = lookupFloat[0]->getLabel();
1059 //get next line to process
1060 //prevent memory leak
1061 for (int i = 0; i < lookupFloat.size(); i++) { delete lookupFloat[i]; }
1062 lookupFloat = input->getSharedRAbundFloatVectors();
1066 if (m->control_pressed) { delete input; return 0; }
1068 //output error messages about any remaining user labels
1069 set<string>::iterator it;
1070 bool needToRun = false;
1071 for (it = userLabels.begin(); it != userLabels.end(); it++) {
1072 m->mothurOut("Your file does not include the label " + *it);
1073 if (processedLabels.count(lastLabel) != 1) {
1074 m->mothurOut(". I will use " + lastLabel + "."); m->mothurOutEndLine();
1077 m->mothurOut(". Please refer to " + lastLabel + "."); m->mothurOutEndLine();
1081 //run last label if you need to
1082 if (needToRun == true) {
1083 for (int i = 0; i < lookupFloat.size(); i++) { if (lookupFloat[i] != NULL) { delete lookupFloat[i]; } }
1084 lookupFloat = input->getSharedRAbundFloatVectors(lastLabel);
1090 catch(exception& e) {
1091 m->errorOut(e, "IndicatorCommand", "getShared");
1095 //**********************************************************************************************************************
1096 vector<float> IndicatorCommand::driver(vector< vector<SharedRAbundFloatVector*> >& groupings, map< vector<int>, vector<int> > groupingsMap, int num, vector<float> indicatorValues, int numIters){
1098 vector<float> pvalues;
1099 pvalues.resize(indicatorValues.size(), 0);
1101 for(int i=0;i<numIters;i++){
1102 if (m->control_pressed) { break; }
1103 groupingsMap = randomizeGroupings(groupings, num);
1104 vector<float> randomIndicatorValues = getValues(groupings, groupingsMap);
1106 for (int j = 0; j < indicatorValues.size(); j++) {
1107 if (randomIndicatorValues[j] >= indicatorValues[j]) { pvalues[j]++; }
1113 }catch(exception& e) {
1114 m->errorOut(e, "IndicatorCommand", "driver");
1118 //**********************************************************************************************************************
1119 vector<float> IndicatorCommand::getPValues(vector< vector<SharedRAbundFloatVector*> >& groupings, map< vector<int>, vector<int> > groupingsMap, int num, vector<float> indicatorValues){
1121 vector<float> pvalues;
1123 #if defined (__APPLE__) || (__MACH__) || (linux) || (__linux)
1124 if(processors == 1){
1125 pvalues = driver(groupings, groupingsMap, num, indicatorValues, iters);
1126 for (int i = 0; i < pvalues.size(); i++) { pvalues[i] /= (double)iters; }
1129 //divide iters between processors
1130 int numItersPerProcessor = iters / processors;
1132 vector<int> processIDS;
1136 //loop through and create all the processes you want
1137 while (process != processors) {
1141 processIDS.push_back(pid); //create map from line number to pid so you can append files in correct order later
1143 }else if (pid == 0){
1144 pvalues = driver(groupings, groupingsMap, num, indicatorValues, numItersPerProcessor);
1146 //pass pvalues to parent
1148 string tempFile = toString(getpid()) + ".pvalues.temp";
1149 m->openOutputFile(tempFile, out);
1152 for (int i = 0; i < pvalues.size(); i++) {
1153 out << pvalues[i] << '\t';
1161 m->mothurOut("[ERROR]: unable to spawn the necessary processes."); m->mothurOutEndLine();
1162 for (int i = 0; i < processIDS.size(); i++) { kill (processIDS[i], SIGINT); }
1168 //special case for last processor in case it doesn't divide evenly
1169 numItersPerProcessor = iters - ((processors-1) * numItersPerProcessor);
1170 pvalues = driver(groupings, groupingsMap, num, indicatorValues, numItersPerProcessor);
1172 //force parent to wait until all the processes are done
1173 for (int i=0;i<processIDS.size();i++) {
1174 int temp = processIDS[i];
1179 for (int i = 0; i < processIDS.size(); i++) {
1181 string tempFile = toString(processIDS[i]) + ".pvalues.temp";
1182 m->openInputFile(tempFile, in);
1184 ////// to do ///////////
1185 int numTemp; numTemp = 0;
1186 for (int j = 0; j < pvalues.size(); j++) {
1187 in >> numTemp; m->gobble(in);
1188 pvalues[j] += numTemp;
1190 in.close(); m->mothurRemove(tempFile);
1192 for (int i = 0; i < pvalues.size(); i++) { pvalues[i] /= (double)iters; }
1195 pvalues = driver(groupings, groupingsMap, num, indicatorValues, iters);
1196 for (int i = 0; i < pvalues.size(); i++) { pvalues[i] /= (double)iters; }
1201 catch(exception& e) {
1202 m->errorOut(e, "IndicatorCommand", "getPValues");
1207 //**********************************************************************************************************************
1208 //same as above, just data type difference
1209 vector<float> IndicatorCommand::driver(vector< vector<SharedRAbundVector*> >& groupings, map< vector<int>, vector<int> > groupingsMap, int num, vector<float> indicatorValues, int numIters){
1211 vector<float> pvalues;
1212 pvalues.resize(indicatorValues.size(), 0);
1214 for(int i=0;i<numIters;i++){
1215 if (m->control_pressed) { break; }
1216 groupingsMap = randomizeGroupings(groupings, num);
1217 vector<float> randomIndicatorValues = getValues(groupings, groupingsMap);
1219 for (int j = 0; j < indicatorValues.size(); j++) {
1220 if (randomIndicatorValues[j] >= indicatorValues[j]) { pvalues[j]++; }
1226 }catch(exception& e) {
1227 m->errorOut(e, "IndicatorCommand", "driver");
1231 //**********************************************************************************************************************
1232 //same as above, just data type difference
1233 vector<float> IndicatorCommand::getPValues(vector< vector<SharedRAbundVector*> >& groupings, map< vector<int>, vector<int> > groupingsMap, int num, vector<float> indicatorValues){
1235 vector<float> pvalues;
1237 #if defined (__APPLE__) || (__MACH__) || (linux) || (__linux)
1238 if(processors == 1){
1239 pvalues = driver(groupings, groupingsMap, num, indicatorValues, iters);
1240 for (int i = 0; i < pvalues.size(); i++) { pvalues[i] /= (double)iters; }
1243 //divide iters between processors
1244 int numItersPerProcessor = iters / processors;
1246 vector<int> processIDS;
1249 //loop through and create all the processes you want
1250 while (process != processors) {
1254 processIDS.push_back(pid); //create map from line number to pid so you can append files in correct order later
1256 }else if (pid == 0){
1257 pvalues = driver(groupings, groupingsMap, num, indicatorValues, numItersPerProcessor);
1259 //pass pvalues to parent
1261 string tempFile = toString(getpid()) + ".pvalues.temp";
1262 m->openOutputFile(tempFile, out);
1265 for (int i = 0; i < pvalues.size(); i++) {
1266 out << pvalues[i] << '\t';
1274 m->mothurOut("[ERROR]: unable to spawn the necessary processes."); m->mothurOutEndLine();
1275 for (int i = 0; i < processIDS.size(); i++) { kill (processIDS[i], SIGINT); }
1281 //special case for last processor in case it doesn't divide evenly
1282 numItersPerProcessor = iters - ((processors-1) * numItersPerProcessor);
1283 pvalues = driver(groupings, groupingsMap, num, indicatorValues, numItersPerProcessor);
1285 //force parent to wait until all the processes are done
1286 for (int i=0;i<processIDS.size();i++) {
1287 int temp = processIDS[i];
1292 for (int i = 0; i < processIDS.size(); i++) {
1294 string tempFile = toString(processIDS[i]) + ".pvalues.temp";
1295 m->openInputFile(tempFile, in);
1297 ////// to do ///////////
1298 int numTemp; numTemp = 0;
1299 for (int j = 0; j < pvalues.size(); j++) {
1300 in >> numTemp; m->gobble(in);
1301 pvalues[j] += numTemp;
1303 in.close(); m->mothurRemove(tempFile);
1305 for (int i = 0; i < pvalues.size(); i++) { pvalues[i] /= (double)iters; }
1308 pvalues = driver(groupings, groupingsMap, num, indicatorValues, iters);
1309 for (int i = 0; i < pvalues.size(); i++) { pvalues[i] /= (double)iters; }
1314 catch(exception& e) {
1315 m->errorOut(e, "IndicatorCommand", "getPValues");
1319 //**********************************************************************************************************************
1320 //swap groups between groupings, in essence randomizing the second column of the design file
1321 map< vector<int>, vector<int> > IndicatorCommand::randomizeGroupings(vector< vector<SharedRAbundVector*> >& groupings, int numLookupGroups){
1324 map< vector<int>, vector<int> > randomGroupings;
1326 for (int i = 0; i < numLookupGroups; i++) {
1327 if (m->control_pressed) {break;}
1329 //get random groups to swap to switch with
1330 //generate random int between 0 and groupings.size()-1
1331 int z = m->getRandomIndex(groupings.size()-1);
1332 int x = m->getRandomIndex(groupings.size()-1);
1333 int a = m->getRandomIndex(groupings[z].size()-1);
1334 int b = m->getRandomIndex(groupings[x].size()-1);
1335 //cout << i << '\t' << z << '\t' << x << '\t' << a << '\t' << b << endl;
1336 //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; }
1341 from.push_back(z); from.push_back(a);
1342 to.push_back(x); to.push_back(b);
1344 randomGroupings[from] = to;
1346 //cout << "done" << endl;
1347 return randomGroupings;
1349 catch(exception& e) {
1350 m->errorOut(e, "IndicatorCommand", "randomizeGroupings");
1354 //**********************************************************************************************************************
1355 //swap groups between groupings, in essence randomizing the second column of the design file
1356 map< vector<int>, vector<int> > IndicatorCommand::randomizeGroupings(vector< vector<SharedRAbundFloatVector*> >& groupings, int numLookupGroups){
1359 map< vector<int>, vector<int> > randomGroupings;
1361 for (int i = 0; i < numLookupGroups; i++) {
1363 //get random groups to swap to switch with
1364 //generate random int between 0 and groupings.size()-1
1365 int z = m->getRandomIndex(groupings.size()-1);
1366 int x = m->getRandomIndex(groupings.size()-1);
1367 int a = m->getRandomIndex(groupings[z].size()-1);
1368 int b = m->getRandomIndex(groupings[x].size()-1);
1369 //cout << i << '\t' << z << '\t' << x << '\t' << a << '\t' << b << endl;
1374 from.push_back(z); from.push_back(a);
1375 to.push_back(x); to.push_back(b);
1377 randomGroupings[from] = to;
1380 return randomGroupings;
1382 catch(exception& e) {
1383 m->errorOut(e, "IndicatorCommand", "randomizeGroupings");
1388 /*****************************************************************/