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();
103 vector<string> tempOutNames;
104 outputTypes["tree"] = tempOutNames;
105 outputTypes["summary"] = tempOutNames;
107 //if the user changes the input directory command factory will send this info to us in the output parameter
108 string inputDir = validParameter.validFile(parameters, "inputdir", false);
109 if (inputDir == "not found"){ inputDir = ""; }
112 it = parameters.find("tree");
113 //user has given a template file
114 if(it != parameters.end()){
115 path = m->hasPath(it->second);
116 //if the user has not given a path then, add inputdir. else leave path alone.
117 if (path == "") { parameters["tree"] = inputDir + it->second; }
120 it = parameters.find("shared");
121 //user has given a template file
122 if(it != parameters.end()){
123 path = m->hasPath(it->second);
124 //if the user has not given a path then, add inputdir. else leave path alone.
125 if (path == "") { parameters["shared"] = inputDir + it->second; }
128 it = parameters.find("relabund");
129 //user has given a template file
130 if(it != parameters.end()){
131 path = m->hasPath(it->second);
132 //if the user has not given a path then, add inputdir. else leave path alone.
133 if (path == "") { parameters["relabund"] = inputDir + it->second; }
136 it = parameters.find("design");
137 //user has given a template file
138 if(it != parameters.end()){
139 path = m->hasPath(it->second);
140 //if the user has not given a path then, add inputdir. else leave path alone.
141 if (path == "") { parameters["design"] = inputDir + it->second; }
145 outputDir = validParameter.validFile(parameters, "outputdir", false); if (outputDir == "not found"){ outputDir = ""; }
147 //check for required parameters
148 treefile = validParameter.validFile(parameters, "tree", true);
149 if (treefile == "not open") { treefile = ""; abort = true; }
150 else if (treefile == "not found") { treefile = ""; }
151 else { m->setTreeFile(treefile); }
153 sharedfile = validParameter.validFile(parameters, "shared", true);
154 if (sharedfile == "not open") { abort = true; }
155 else if (sharedfile == "not found") { sharedfile = ""; }
156 else { inputFileName = sharedfile; m->setSharedFile(sharedfile); }
158 relabundfile = validParameter.validFile(parameters, "relabund", true);
159 if (relabundfile == "not open") { abort = true; }
160 else if (relabundfile == "not found") { relabundfile = ""; }
161 else { inputFileName = relabundfile; m->setRelAbundFile(relabundfile); }
163 designfile = validParameter.validFile(parameters, "design", true);
164 if (designfile == "not open") { designfile = ""; abort = true; }
165 else if (designfile == "not found") { designfile = ""; }
166 else { m->setDesignFile(designfile); }
168 groups = validParameter.validFile(parameters, "groups", false);
169 if (groups == "not found") { groups = ""; Groups.push_back("all"); }
170 else { m->splitAtDash(groups, Groups); }
171 m->setGroups(Groups);
173 label = validParameter.validFile(parameters, "label", false);
174 if (label == "not found") { label = ""; m->mothurOut("You did not provide a label, I will use the first label in your inputfile."); m->mothurOutEndLine(); label=""; }
176 string temp = validParameter.validFile(parameters, "iters", false); if (temp == "not found") { temp = "1000"; }
177 m->mothurConvert(temp, iters);
179 temp = validParameter.validFile(parameters, "processors", false); if (temp == "not found"){ temp = m->getProcessors(); }
180 m->setProcessors(temp);
181 m->mothurConvert(temp, processors);
183 if ((relabundfile == "") && (sharedfile == "")) {
184 //is there are current file available for either of these?
185 //give priority to shared, then relabund
186 sharedfile = m->getSharedFile();
187 if (sharedfile != "") { inputFileName = sharedfile; m->mothurOut("Using " + sharedfile + " as input file for the shared parameter."); m->mothurOutEndLine(); }
189 relabundfile = m->getRelAbundFile();
190 if (relabundfile != "") { inputFileName = relabundfile; m->mothurOut("Using " + relabundfile + " as input file for the relabund parameter."); m->mothurOutEndLine(); }
192 m->mothurOut("No valid current files. You must provide a shared or relabund."); m->mothurOutEndLine();
199 if ((designfile == "") && (treefile == "")) {
200 treefile = m->getTreeFile();
201 if (treefile != "") { m->mothurOut("Using " + treefile + " as input file for the tree parameter."); m->mothurOutEndLine(); }
203 designfile = m->getDesignFile();
204 if (designfile != "") { m->mothurOut("Using " + designfile + " as input file for the design parameter."); m->mothurOutEndLine(); }
206 m->mothurOut("[ERROR]: You must provide either a tree or design file."); m->mothurOutEndLine(); abort = true;
211 if ((relabundfile != "") && (sharedfile != "")) { m->mothurOut("[ERROR]: You may not use both a shared and relabund file."); m->mothurOutEndLine(); abort = true; }
216 catch(exception& e) {
217 m->errorOut(e, "IndicatorCommand", "IndicatorCommand");
221 //**********************************************************************************************************************
223 int IndicatorCommand::execute(){
226 if (abort == true) { if (calledHelp) { return 0; } return 2; }
228 cout.setf(ios::fixed, ios::floatfield); cout.setf(ios::showpoint);
230 int start = time(NULL);
232 //read designfile if given and set up groups for read of sharedfiles
233 if (designfile != "") {
234 designMap = new GroupMap(designfile);
235 designMap->readDesignMap();
237 //fill Groups - checks for "all" and for any typo groups
239 vector<string> nameGroups = designMap->getNamesOfGroups();
240 util.setGroups(Groups, nameGroups);
241 designMap->setNamesOfGroups(nameGroups);
243 vector<string> namesSeqs = designMap->getNamesSeqs(Groups);
244 m->setGroups(namesSeqs);
247 /***************************************************/
248 // use smart distancing to get right sharedRabund //
249 /***************************************************/
250 if (sharedfile != "") {
252 if (m->control_pressed) { if (designfile != "") { delete designMap; } for (int i = 0; i < lookup.size(); i++) { delete lookup[i]; } return 0; }
253 if (lookup[0] == NULL) { m->mothurOut("[ERROR] reading shared file."); m->mothurOutEndLine(); return 0; }
256 if (m->control_pressed) { if (designfile != "") { delete designMap; } for (int i = 0; i < lookupFloat.size(); i++) { delete lookupFloat[i]; } return 0; }
257 if (lookupFloat[0] == NULL) { m->mothurOut("[ERROR] reading relabund file."); m->mothurOutEndLine(); return 0; }
260 //reset groups if needed
261 if (designfile != "") { m->setGroups(Groups); }
263 /***************************************************/
264 // reading tree info //
265 /***************************************************/
266 if (treefile != "") {
267 string groupfile = "";
268 m->setTreeFile(treefile);
269 Tree* tree = new Tree(treefile); delete tree; //extracts names from tree to make faked out groupmap
270 treeMap = new TreeMap();
271 bool mismatch = false;
273 for (int i = 0; i < m->Treenames.size(); i++) {
274 //sanity check - is this a group that is not in the sharedfile?
275 if (designfile == "") {
276 if (!(m->inUsersGroups(m->Treenames[i], m->getAllGroups()))) {
277 m->mothurOut("[ERROR]: " + m->Treenames[i] + " is not a group in your shared or relabund file."); m->mothurOutEndLine();
280 treeMap->addSeq(m->Treenames[i], "Group1");
282 vector<string> myGroups; myGroups.push_back(m->Treenames[i]);
283 vector<string> myNames = designMap->getNamesSeqs(myGroups);
285 for(int k = 0; k < myNames.size(); k++) {
286 if (!(m->inUsersGroups(myNames[k], m->getAllGroups()))) {
287 m->mothurOut("[ERROR]: " + myNames[k] + " is not a group in your shared or relabund file."); m->mothurOutEndLine();
291 treeMap->addSeq(m->Treenames[i], "Group1");
295 if ((designfile != "") && (m->Treenames.size() != Groups.size())) { cout << Groups.size() << '\t' << m->Treenames.size() << endl; m->mothurOut("[ERROR]: You design file does not match your tree, aborting."); m->mothurOutEndLine(); mismatch = true; }
297 if (mismatch) { //cleanup and exit
298 if (designfile != "") { delete designMap; }
299 if (sharedfile != "") { for (int i = 0; i < lookup.size(); i++) { delete lookup[i]; } }
300 else { for (int i = 0; i < lookupFloat.size(); i++) { delete lookupFloat[i]; } }
305 read = new ReadNewickTree(treefile);
306 int readOk = read->read(treeMap);
308 if (readOk != 0) { m->mothurOut("Read Terminated."); m->mothurOutEndLine(); delete treeMap; delete read; return 0; }
310 vector<Tree*> T = read->getTrees();
314 if (m->control_pressed) {
315 if (designfile != "") { delete designMap; }
316 if (sharedfile != "") { for (int i = 0; i < lookup.size(); i++) { delete lookup[i]; } }
317 else { for (int i = 0; i < lookupFloat.size(); i++) { delete lookupFloat[i]; } }
318 for (int i = 0; i < T.size(); i++) { delete T[i]; } delete treeMap; return 0;
321 map<string, string> nameMap;
322 T[0]->assembleTree(nameMap);
324 /***************************************************/
325 // create ouptut tree - respecting pickedGroups //
326 /***************************************************/
327 Tree* outputTree = new Tree(m->getNumGroups(), treeMap);
329 outputTree->getSubTree(T[0], m->getGroups());
330 outputTree->assembleTree(nameMap);
332 //no longer need original tree, we have output tree to use and label
333 for (int i = 0; i < T.size(); i++) { delete T[i]; }
335 if (m->control_pressed) {
336 if (designfile != "") { delete designMap; }
337 if (sharedfile != "") { for (int i = 0; i < lookup.size(); i++) { delete lookup[i]; } }
338 else { for (int i = 0; i < lookupFloat.size(); i++) { delete lookupFloat[i]; } }
339 delete outputTree; delete treeMap; return 0;
342 /***************************************************/
343 // get indicator species values //
344 /***************************************************/
345 GetIndicatorSpecies(outputTree);
346 delete outputTree; delete treeMap;
348 }else { //run with design file only
349 //get indicator species
350 GetIndicatorSpecies();
353 if (designfile != "") { delete designMap; }
354 if (sharedfile != "") { for (int i = 0; i < lookup.size(); i++) { delete lookup[i]; } }
355 else { for (int i = 0; i < lookupFloat.size(); i++) { delete lookupFloat[i]; } }
357 if (m->control_pressed) { for (int i = 0; i < outputNames.size(); i++) { m->mothurRemove(outputNames[i]); } return 0; }
359 //set tree file as new current treefile
360 if (treefile != "") {
362 itTypes = outputTypes.find("tree");
363 if (itTypes != outputTypes.end()) {
364 if ((itTypes->second).size() != 0) { current = (itTypes->second)[0]; m->setTreeFile(current); }
368 m->mothurOutEndLine(); m->mothurOutEndLine(); m->mothurOut("It took " + toString(time(NULL) - start) + " secs to find the indicator species."); m->mothurOutEndLine();
369 m->mothurOutEndLine();
370 m->mothurOut("Output File Names: "); m->mothurOutEndLine();
371 for (int i = 0; i < outputNames.size(); i++) { m->mothurOut(outputNames[i]); m->mothurOutEndLine(); }
372 m->mothurOutEndLine();
376 catch(exception& e) {
377 m->errorOut(e, "IndicatorCommand", "execute");
381 //**********************************************************************************************************************
382 //divide shared or relabund file by groupings in the design file
383 //report all otu values to file
384 int IndicatorCommand::GetIndicatorSpecies(){
386 string thisOutputDir = outputDir;
387 if (outputDir == "") { thisOutputDir += m->hasPath(inputFileName); }
388 string outputFileName = thisOutputDir + m->getRootName(m->getSimpleName(inputFileName)) + "indicator.summary";
389 outputNames.push_back(outputFileName); outputTypes["summary"].push_back(outputFileName);
392 m->openOutputFile(outputFileName, out);
393 out.setf(ios::fixed, ios::floatfield); out.setf(ios::showpoint);
394 m->mothurOutEndLine(); m->mothurOut("Species\tIndicatorValue\tpValue\n");
397 if (sharedfile != "") { numBins = lookup[0]->getNumBins(); }
398 else { numBins = lookupFloat[0]->getNumBins(); }
400 if (m->control_pressed) { out.close(); return 0; }
402 /*****************************************************/
403 //create vectors containing rabund info //
404 /*****************************************************/
406 vector<float> indicatorValues; //size of numBins
407 vector<float> pValues;
408 map< vector<int>, vector<int> > randomGroupingsMap; //maps location in groupings to location in groupings, ie, [0][0] -> [1][2]. This is so we don't have to actually move the sharedRabundVectors.
410 if (sharedfile != "") {
411 vector< vector<SharedRAbundVector*> > groupings;
412 set<string> groupsAlreadyAdded;
413 vector<SharedRAbundVector*> subset;
416 for (int i = 0; i < (designMap->getNamesOfGroups()).size(); i++) {
418 for (int k = 0; k < lookup.size(); k++) {
419 //are you from this grouping?
420 if (designMap->getGroup(lookup[k]->getGroup()) == (designMap->getNamesOfGroups())[i]) {
421 subset.push_back(lookup[k]);
422 groupsAlreadyAdded.insert(lookup[k]->getGroup());
425 if (subset.size() != 0) { groupings.push_back(subset); }
429 if (groupsAlreadyAdded.size() != lookup.size()) { m->mothurOut("[ERROR]: could not make proper groupings."); m->mothurOutEndLine(); }
431 indicatorValues = getValues(groupings, randomGroupingsMap);
433 pValues = getPValues(groupings, randomGroupingsMap, lookup.size(), indicatorValues);
435 vector< vector<SharedRAbundFloatVector*> > groupings;
436 set<string> groupsAlreadyAdded;
437 vector<SharedRAbundFloatVector*> subset;
440 for (int i = 0; i < (designMap->getNamesOfGroups()).size(); i++) {
441 for (int k = 0; k < lookupFloat.size(); k++) {
442 //are you from this grouping?
443 if (designMap->getGroup(lookupFloat[k]->getGroup()) == (designMap->getNamesOfGroups())[i]) {
444 subset.push_back(lookupFloat[k]);
445 groupsAlreadyAdded.insert(lookupFloat[k]->getGroup());
448 if (subset.size() != 0) { groupings.push_back(subset); }
452 if (groupsAlreadyAdded.size() != lookupFloat.size()) { m->mothurOut("[ERROR]: could not make proper groupings."); m->mothurOutEndLine(); }
454 indicatorValues = getValues(groupings, randomGroupingsMap);
456 pValues = getPValues(groupings, randomGroupingsMap, lookupFloat.size(), indicatorValues);
459 if (m->control_pressed) { out.close(); return 0; }
462 /******************************************************/
463 //output indicator values to table form //
464 /*****************************************************/
465 out << "OTU\tIndicator_Value\tpValue" << endl;
466 for (int j = 0; j < indicatorValues.size(); j++) {
468 if (m->control_pressed) { out.close(); return 0; }
470 out << m->currentBinLabels[j] << '\t' << indicatorValues[j] << '\t';
472 if (pValues[j] > (1/(float)iters)) { out << pValues[j] << endl; }
473 else { out << "<" << (1/(float)iters) << endl; }
475 if (pValues[j] <= 0.05) {
476 cout << m->currentBinLabels[j] << '\t' << indicatorValues[j] << '\t';
477 string pValueString = "<" + toString((1/(float)iters));
478 if (pValues[j] > (1/(float)iters)) { pValueString = toString(pValues[j]); cout << pValues[j];}
479 else { cout << "<" << (1/(float)iters); }
480 m->mothurOutJustToLog(m->currentBinLabels[j] + "\t" + toString(indicatorValues[j]) + "\t" + pValueString);
481 m->mothurOutEndLine();
489 catch(exception& e) {
490 m->errorOut(e, "IndicatorCommand", "GetIndicatorSpecies");
494 //**********************************************************************************************************************
495 //traverse tree finding indicator species values for each otu at each node
496 //label node with otu number that has highest indicator value
497 //report all otu values to file
498 int IndicatorCommand::GetIndicatorSpecies(Tree*& T){
501 string thisOutputDir = outputDir;
502 if (outputDir == "") { thisOutputDir += m->hasPath(inputFileName); }
503 string outputFileName = thisOutputDir + m->getRootName(m->getSimpleName(inputFileName)) + "indicator.summary";
504 outputNames.push_back(outputFileName); outputTypes["summary"].push_back(outputFileName);
507 m->openOutputFile(outputFileName, out);
508 out.setf(ios::fixed, ios::floatfield); out.setf(ios::showpoint);
511 if (sharedfile != "") { numBins = lookup[0]->getNumBins(); }
512 else { numBins = lookupFloat[0]->getNumBins(); }
516 for (int i = 0; i < numBins; i++) { out << m->currentBinLabels[i] << "_IndValue" << '\t' << "pValue" << '\t'; }
519 m->mothurOutEndLine(); m->mothurOut("Node\tSpecies\tIndicatorValue\tpValue\n");
521 string treeOutputDir = outputDir;
522 if (outputDir == "") { treeOutputDir += m->hasPath(treefile); }
523 string outputTreeFileName = treeOutputDir + m->getRootName(m->getSimpleName(treefile)) + "indicator.tre";
526 //create a map from tree node index to names of descendants, save time later to know which sharedRabund you need
527 map<int, set<string> > nodeToDescendants;
528 map<int, set<int> > descendantNodes;
529 for (int i = 0; i < T->getNumNodes(); i++) {
530 if (m->control_pressed) { return 0; }
532 nodeToDescendants[i] = getDescendantList(T, i, nodeToDescendants, descendantNodes);
535 //you need the distances to leaf to decide grouping below
536 //this will also set branch lengths if the tree does not include them
537 map<int, float> distToRoot = getDistToRoot(T);
540 for (int i = T->getNumLeaves(); i < T->getNumNodes(); i++) {
541 //cout << endl << i+1 << endl;
542 if (m->control_pressed) { out.close(); return 0; }
544 /*****************************************************/
545 //create vectors containing rabund info //
546 /*****************************************************/
548 vector<float> indicatorValues; //size of numBins
549 vector<float> pValues;
550 map< vector<int>, vector<int> > randomGroupingsMap; //maps location in groupings to location in groupings, ie, [0][0] -> [1][2]. This is so we don't have to actually move the sharedRabundVectors.
552 if (sharedfile != "") {
553 vector< vector<SharedRAbundVector*> > groupings;
555 //get nodes that will be a valid grouping
556 //you are valid if you are not one of my descendants
557 //AND your distToRoot is >= mine
558 //AND you were not added as part of a larger grouping. Largest nodes are added first.
560 set<string> groupsAlreadyAdded;
561 //create a grouping with my grouping
562 vector<SharedRAbundVector*> subset;
564 int doneCount = nodeToDescendants[i].size();
565 for (int k = 0; k < lookup.size(); k++) {
566 //is this descendant of i
567 if ((nodeToDescendants[i].count(lookup[k]->getGroup()) != 0)) {
568 subset.push_back(lookup[k]);
569 groupsAlreadyAdded.insert(lookup[k]->getGroup());
572 if (count == doneCount) { break; } //quit once you get the rabunds for this grouping
574 if (subset.size() != 0) { groupings.push_back(subset); }
577 for (int j = (T->getNumNodes()-1); j >= 0; j--) {
580 if ((descendantNodes[i].count(j) == 0) && (distToRoot[j] >= distToRoot[i])) {
581 vector<SharedRAbundVector*> subset;
583 int doneCount = nodeToDescendants[j].size();
584 for (int k = 0; k < lookup.size(); k++) {
585 //is this descendant of j, and we didn't already add this as part of a larger grouping
586 if ((nodeToDescendants[j].count(lookup[k]->getGroup()) != 0) && (groupsAlreadyAdded.count(lookup[k]->getGroup()) == 0)) {
587 subset.push_back(lookup[k]);
588 groupsAlreadyAdded.insert(lookup[k]->getGroup());
591 if (count == doneCount) { break; } //quit once you get the rabunds for this grouping
594 //if subset.size == 0 then the node was added as part of a larger grouping
595 if (subset.size() != 0) { groupings.push_back(subset); }
599 if (groupsAlreadyAdded.size() != lookup.size()) { m->mothurOut("[ERROR]: could not make proper groupings."); m->mothurOutEndLine(); }
601 indicatorValues = getValues(groupings, randomGroupingsMap);
603 pValues = getPValues(groupings, randomGroupingsMap, lookup.size(), indicatorValues);
605 vector< vector<SharedRAbundFloatVector*> > groupings;
607 //get nodes that will be a valid grouping
608 //you are valid if you are not one of my descendants
609 //AND your distToRoot is >= mine
610 //AND you were not added as part of a larger grouping. Largest nodes are added first.
612 set<string> groupsAlreadyAdded;
613 //create a grouping with my grouping
614 vector<SharedRAbundFloatVector*> subset;
616 int doneCount = nodeToDescendants[i].size();
617 for (int k = 0; k < lookupFloat.size(); k++) {
618 //is this descendant of i
619 if ((nodeToDescendants[i].count(lookupFloat[k]->getGroup()) != 0)) {
620 subset.push_back(lookupFloat[k]);
621 groupsAlreadyAdded.insert(lookupFloat[k]->getGroup());
624 if (count == doneCount) { break; } //quit once you get the rabunds for this grouping
626 if (subset.size() != 0) { groupings.push_back(subset); }
628 for (int j = (T->getNumNodes()-1); j >= 0; j--) {
629 if ((descendantNodes[i].count(j) == 0) && (distToRoot[j] >= distToRoot[i])) {
630 vector<SharedRAbundFloatVector*> subset;
632 int doneCount = nodeToDescendants[j].size();
633 for (int k = 0; k < lookupFloat.size(); k++) {
634 //is this descendant of j, and we didn't already add this as part of a larger grouping
635 if ((nodeToDescendants[j].count(lookupFloat[k]->getGroup()) != 0) && (groupsAlreadyAdded.count(lookupFloat[k]->getGroup()) == 0)) {
636 subset.push_back(lookupFloat[k]);
637 groupsAlreadyAdded.insert(lookupFloat[k]->getGroup());
640 if (count == doneCount) { break; } //quit once you get the rabunds for this grouping
643 //if subset.size == 0 then the node was added as part of a larger grouping
644 if (subset.size() != 0) { groupings.push_back(subset); }
648 if (groupsAlreadyAdded.size() != lookupFloat.size()) { m->mothurOut("[ERROR]: could not make proper groupings."); m->mothurOutEndLine(); }
650 indicatorValues = getValues(groupings, randomGroupingsMap);
652 pValues = getPValues(groupings, randomGroupingsMap, lookupFloat.size(), indicatorValues);
655 if (m->control_pressed) { out.close(); return 0; }
658 /******************************************************/
659 //output indicator values to table form + label tree //
660 /*****************************************************/
661 out << (i+1) << '\t';
662 for (int j = 0; j < indicatorValues.size(); j++) {
664 if (m->control_pressed) { out.close(); return 0; }
666 if (pValues[j] < (1/(float)iters)) {
667 out << indicatorValues[j] << '\t' << '<' << (1/(float)iters) << '\t';
669 out << indicatorValues[j] << '\t' << pValues[j] << '\t';
672 if (pValues[j] <= 0.05) {
673 cout << i+1 << '\t' << m->currentBinLabels[j] << '\t' << indicatorValues[j] << '\t';
674 string pValueString = "<" + toString((1/(float)iters));
675 if (pValues[j] > (1/(float)iters)) { pValueString = toString(pValues[j]); cout << pValues[j];}
676 else { cout << "<" << (1/(float)iters); }
677 m->mothurOutJustToLog(toString(i) + "\t" + m->currentBinLabels[j] + "\t" + toString(indicatorValues[j]) + "\t" + pValueString);
678 m->mothurOutEndLine();
683 T->tree[i].setLabel((i+1));
688 m->openOutputFile(outputTreeFileName, outTree);
689 outputNames.push_back(outputTreeFileName); outputTypes["tree"].push_back(outputTreeFileName);
691 T->print(outTree, "both");
696 catch(exception& e) {
697 m->errorOut(e, "IndicatorCommand", "GetIndicatorSpecies");
701 //**********************************************************************************************************************
702 vector<float> IndicatorCommand::getValues(vector< vector<SharedRAbundFloatVector*> >& groupings, map< vector<int>, vector<int> > groupingsMap){
704 vector<float> values;
705 map< vector<int>, vector<int> >::iterator it;
708 for (int i = 0; i < groupings[0][0]->getNumBins(); i++) {
710 if (m->control_pressed) { return values; }
713 float AijDenominator = 0.0;
716 //get overall abundance of each grouping
717 for (int j = 0; j < groupings.size(); j++) {
719 float totalAbund = 0;
721 for (int k = 0; k < groupings[j].size(); k++) {
722 vector<int> temp; temp.push_back(j); temp.push_back(k);
723 it = groupingsMap.find(temp);
725 if (it == groupingsMap.end()) { //this one didnt get moved
726 totalAbund += groupings[j][k]->getAbundance(i);
727 if (groupings[j][k]->getAbundance(i) != 0.0) { numNotZero++; }
729 totalAbund += groupings[(it->second)[0]][(it->second)[1]]->getAbundance(i);
730 if (groupings[(it->second)[0]][(it->second)[1]]->getAbundance(i) != 0.0) { numNotZero++; }
736 float Aij = (totalAbund / (float) groupings[j].size());
737 terms.push_back(Aij);
739 //percentage of sites represented
740 Bij.push_back(numNotZero / (float) groupings[j].size());
742 AijDenominator += Aij;
745 float maxIndVal = 0.0;
746 for (int j = 0; j < terms.size(); j++) {
747 float thisAij = (terms[j] / AijDenominator); //relative abundance
748 float thisValue = thisAij * Bij[j] * 100.0;
751 if (thisValue > maxIndVal) { maxIndVal = thisValue; }
754 values.push_back(maxIndVal);
759 catch(exception& e) {
760 m->errorOut(e, "IndicatorCommand", "getValues");
764 //**********************************************************************************************************************
765 //same as above, just data type difference
766 vector<float> IndicatorCommand::getValues(vector< vector<SharedRAbundVector*> >& groupings, map< vector<int>, vector<int> > groupingsMap){
768 vector<float> values;
770 /*for (int j = 0; j < groupings.size(); j++) {
771 cout << "grouping " << j << endl;
772 for (int k = 0; k < groupings[j].size(); k++) {
773 cout << groupings[j][k]->getGroup() << endl;
776 map< vector<int>, vector<int> >::iterator it;
779 for (int i = 0; i < groupings[0][0]->getNumBins(); i++) {
781 float AijDenominator = 0.0;
783 //get overall abundance of each grouping
784 for (int j = 0; j < groupings.size(); j++) {
786 int totalAbund = 0.0;
788 for (int k = 0; k < groupings[j].size(); k++) {
789 vector<int> temp; temp.push_back(j); temp.push_back(k);
790 it = groupingsMap.find(temp);
792 if (it == groupingsMap.end()) { //this one didnt get moved
793 totalAbund += groupings[j][k]->getAbundance(i);
794 if (groupings[j][k]->getAbundance(i) != 0.0) { numNotZero++; }
796 totalAbund += groupings[(it->second)[0]][(it->second)[1]]->getAbundance(i);
797 if (groupings[(it->second)[0]][(it->second)[1]]->getAbundance(i) != 0.0) { numNotZero++; }
803 float Aij = (totalAbund / (float) groupings[j].size());
804 terms.push_back(Aij);
806 //percentage of sites represented
807 Bij.push_back(numNotZero / (float) groupings[j].size());
809 AijDenominator += Aij;
812 float maxIndVal = 0.0;
813 for (int j = 0; j < terms.size(); j++) {
814 float thisAij = (terms[j] / AijDenominator); //relative abundance
815 float thisValue = thisAij * Bij[j] * 100.0;
818 if (thisValue > maxIndVal) { maxIndVal = thisValue; }
821 values.push_back(maxIndVal);
826 catch(exception& e) {
827 m->errorOut(e, "IndicatorCommand", "getValues");
831 //**********************************************************************************************************************
832 //you need the distances to root to decide groupings
833 //this will also set branch lengths if the tree does not include them
834 map<int, float> IndicatorCommand::getDistToRoot(Tree*& T){
836 map<int, float> dists;
838 bool hasBranchLengths = false;
839 for (int i = 0; i < T->getNumNodes(); i++) {
840 if (T->tree[i].getBranchLength() > 0.0) { hasBranchLengths = true; break; }
843 //set branchlengths if needed
844 if (!hasBranchLengths) {
845 for (int i = 0; i < T->getNumNodes(); i++) {
847 int lc = T->tree[i].getLChild();
848 int rc = T->tree[i].getRChild();
850 if (lc == -1) { // you are a leaf
851 //if you are a leaf set you priliminary length to 1.0, this may adjust later
852 T->tree[i].setBranchLength(1.0);
854 }else{ // you are an internal node
855 //look at your children's length to leaf
856 float ldist = dists[lc];
857 float rdist = dists[rc];
859 float greater = ldist;
860 if (rdist > greater) { greater = rdist; dists[i] = ldist + 1.0;}
861 else { dists[i] = rdist + 1.0; }
864 //branch length = difference + 1
865 T->tree[lc].setBranchLength((abs(ldist-greater) + 1.0));
866 T->tree[rc].setBranchLength((abs(rdist-greater) + 1.0));
873 for (int i = 0; i < T->getNumNodes(); i++) {
878 while(T->tree[index].getParent() != -1){
879 if (T->tree[index].getBranchLength() != -1) {
880 sum += abs(T->tree[index].getBranchLength());
882 index = T->tree[index].getParent();
890 catch(exception& e) {
891 m->errorOut(e, "IndicatorCommand", "getLengthToLeaf");
895 //**********************************************************************************************************************
896 set<string> IndicatorCommand::getDescendantList(Tree*& T, int i, map<int, set<string> > descendants, map<int, set<int> >& nodes){
900 set<string>::iterator it;
902 int lc = T->tree[i].getLChild();
903 int rc = T->tree[i].getRChild();
905 if (lc == -1) { //you are a leaf your only descendant is yourself
906 set<int> temp; temp.insert(i);
909 if (designfile == "") {
910 names.insert(T->tree[i].getName());
912 vector<string> myGroup; myGroup.push_back(T->tree[i].getName());
913 vector<string> myReps = designMap->getNamesSeqs(myGroup);
914 for (int k = 0; k < myReps.size(); k++) {
915 names.insert(myReps[k]);
919 }else{ //your descedants are the combination of your childrens descendants
920 names = descendants[lc];
921 nodes[i] = nodes[lc];
922 for (it = descendants[rc].begin(); it != descendants[rc].end(); it++) {
925 for (set<int>::iterator itNum = nodes[rc].begin(); itNum != nodes[rc].end(); itNum++) {
926 nodes[i].insert(*itNum);
928 //you are your own descendant
934 catch(exception& e) {
935 m->errorOut(e, "IndicatorCommand", "getDescendantList");
939 //**********************************************************************************************************************
940 int IndicatorCommand::getShared(){
942 InputData* input = new InputData(sharedfile, "sharedfile");
943 lookup = input->getSharedRAbundVectors();
944 string lastLabel = lookup[0]->getLabel();
946 if (label == "") { label = lastLabel; delete input; return 0; }
948 //if the users enters label "0.06" and there is no "0.06" in their file use the next lowest label.
949 set<string> labels; labels.insert(label);
950 set<string> processedLabels;
951 set<string> userLabels = labels;
953 //as long as you are not at the end of the file or done wih the lines you want
954 while((lookup[0] != NULL) && (userLabels.size() != 0)) {
955 if (m->control_pressed) { delete input; return 0; }
957 if(labels.count(lookup[0]->getLabel()) == 1){
958 processedLabels.insert(lookup[0]->getLabel());
959 userLabels.erase(lookup[0]->getLabel());
963 if ((m->anyLabelsToProcess(lookup[0]->getLabel(), userLabels, "") == true) && (processedLabels.count(lastLabel) != 1)) {
964 string saveLabel = lookup[0]->getLabel();
966 for (int i = 0; i < lookup.size(); i++) { delete lookup[i]; }
967 lookup = input->getSharedRAbundVectors(lastLabel);
969 processedLabels.insert(lookup[0]->getLabel());
970 userLabels.erase(lookup[0]->getLabel());
972 //restore real lastlabel to save below
973 lookup[0]->setLabel(saveLabel);
977 lastLabel = lookup[0]->getLabel();
979 //get next line to process
980 //prevent memory leak
981 for (int i = 0; i < lookup.size(); i++) { delete lookup[i]; }
982 lookup = input->getSharedRAbundVectors();
986 if (m->control_pressed) { delete input; return 0; }
988 //output error messages about any remaining user labels
989 set<string>::iterator it;
990 bool needToRun = false;
991 for (it = userLabels.begin(); it != userLabels.end(); it++) {
992 m->mothurOut("Your file does not include the label " + *it);
993 if (processedLabels.count(lastLabel) != 1) {
994 m->mothurOut(". I will use " + lastLabel + "."); m->mothurOutEndLine();
997 m->mothurOut(". Please refer to " + lastLabel + "."); m->mothurOutEndLine();
1001 //run last label if you need to
1002 if (needToRun == true) {
1003 for (int i = 0; i < lookup.size(); i++) { if (lookup[i] != NULL) { delete lookup[i]; } }
1004 lookup = input->getSharedRAbundVectors(lastLabel);
1010 catch(exception& e) {
1011 m->errorOut(e, "IndicatorCommand", "getShared");
1015 //**********************************************************************************************************************
1016 int IndicatorCommand::getSharedFloat(){
1018 InputData* input = new InputData(relabundfile, "relabund");
1019 lookupFloat = input->getSharedRAbundFloatVectors();
1020 string lastLabel = lookupFloat[0]->getLabel();
1022 if (label == "") { label = lastLabel; delete input; return 0; }
1024 //if the users enters label "0.06" and there is no "0.06" in their file use the next lowest label.
1025 set<string> labels; labels.insert(label);
1026 set<string> processedLabels;
1027 set<string> userLabels = labels;
1029 //as long as you are not at the end of the file or done wih the lines you want
1030 while((lookupFloat[0] != NULL) && (userLabels.size() != 0)) {
1032 if (m->control_pressed) { delete input; return 0; }
1034 if(labels.count(lookupFloat[0]->getLabel()) == 1){
1035 processedLabels.insert(lookupFloat[0]->getLabel());
1036 userLabels.erase(lookupFloat[0]->getLabel());
1040 if ((m->anyLabelsToProcess(lookupFloat[0]->getLabel(), userLabels, "") == true) && (processedLabels.count(lastLabel) != 1)) {
1041 string saveLabel = lookupFloat[0]->getLabel();
1043 for (int i = 0; i < lookupFloat.size(); i++) { delete lookupFloat[i]; }
1044 lookupFloat = input->getSharedRAbundFloatVectors(lastLabel);
1046 processedLabels.insert(lookupFloat[0]->getLabel());
1047 userLabels.erase(lookupFloat[0]->getLabel());
1049 //restore real lastlabel to save below
1050 lookupFloat[0]->setLabel(saveLabel);
1054 lastLabel = lookupFloat[0]->getLabel();
1056 //get next line to process
1057 //prevent memory leak
1058 for (int i = 0; i < lookupFloat.size(); i++) { delete lookupFloat[i]; }
1059 lookupFloat = input->getSharedRAbundFloatVectors();
1063 if (m->control_pressed) { delete input; return 0; }
1065 //output error messages about any remaining user labels
1066 set<string>::iterator it;
1067 bool needToRun = false;
1068 for (it = userLabels.begin(); it != userLabels.end(); it++) {
1069 m->mothurOut("Your file does not include the label " + *it);
1070 if (processedLabels.count(lastLabel) != 1) {
1071 m->mothurOut(". I will use " + lastLabel + "."); m->mothurOutEndLine();
1074 m->mothurOut(". Please refer to " + lastLabel + "."); m->mothurOutEndLine();
1078 //run last label if you need to
1079 if (needToRun == true) {
1080 for (int i = 0; i < lookupFloat.size(); i++) { if (lookupFloat[i] != NULL) { delete lookupFloat[i]; } }
1081 lookupFloat = input->getSharedRAbundFloatVectors(lastLabel);
1087 catch(exception& e) {
1088 m->errorOut(e, "IndicatorCommand", "getShared");
1092 //**********************************************************************************************************************
1093 vector<float> IndicatorCommand::driver(vector< vector<SharedRAbundFloatVector*> >& groupings, map< vector<int>, vector<int> > groupingsMap, int num, vector<float> indicatorValues, int numIters){
1095 vector<float> pvalues;
1096 pvalues.resize(indicatorValues.size(), 0);
1098 for(int i=0;i<numIters;i++){
1099 if (m->control_pressed) { break; }
1100 groupingsMap = randomizeGroupings(groupings, num);
1101 vector<float> randomIndicatorValues = getValues(groupings, groupingsMap);
1103 for (int j = 0; j < indicatorValues.size(); j++) {
1104 if (randomIndicatorValues[j] >= indicatorValues[j]) { pvalues[j]++; }
1110 }catch(exception& e) {
1111 m->errorOut(e, "IndicatorCommand", "driver");
1115 //**********************************************************************************************************************
1116 vector<float> IndicatorCommand::getPValues(vector< vector<SharedRAbundFloatVector*> >& groupings, map< vector<int>, vector<int> > groupingsMap, int num, vector<float> indicatorValues){
1118 vector<float> pvalues;
1120 #if defined (__APPLE__) || (__MACH__) || (linux) || (__linux) || (__linux__) || (__unix__) || (__unix)
1121 if(processors == 1){
1122 pvalues = driver(groupings, groupingsMap, num, indicatorValues, iters);
1123 for (int i = 0; i < pvalues.size(); i++) { pvalues[i] /= (double)iters; }
1126 //divide iters between processors
1127 int numItersPerProcessor = iters / processors;
1129 vector<int> processIDS;
1133 //loop through and create all the processes you want
1134 while (process != processors) {
1138 processIDS.push_back(pid); //create map from line number to pid so you can append files in correct order later
1140 }else if (pid == 0){
1141 pvalues = driver(groupings, groupingsMap, num, indicatorValues, numItersPerProcessor);
1143 //pass pvalues to parent
1145 string tempFile = toString(getpid()) + ".pvalues.temp";
1146 m->openOutputFile(tempFile, out);
1149 for (int i = 0; i < pvalues.size(); i++) {
1150 out << pvalues[i] << '\t';
1158 m->mothurOut("[ERROR]: unable to spawn the necessary processes."); m->mothurOutEndLine();
1159 for (int i = 0; i < processIDS.size(); i++) { kill (processIDS[i], SIGINT); }
1165 //special case for last processor in case it doesn't divide evenly
1166 numItersPerProcessor = iters - ((processors-1) * numItersPerProcessor);
1167 pvalues = driver(groupings, groupingsMap, num, indicatorValues, numItersPerProcessor);
1169 //force parent to wait until all the processes are done
1170 for (int i=0;i<processIDS.size();i++) {
1171 int temp = processIDS[i];
1176 for (int i = 0; i < processIDS.size(); i++) {
1178 string tempFile = toString(processIDS[i]) + ".pvalues.temp";
1179 m->openInputFile(tempFile, in);
1181 ////// to do ///////////
1182 int numTemp; numTemp = 0;
1183 for (int j = 0; j < pvalues.size(); j++) {
1184 in >> numTemp; m->gobble(in);
1185 pvalues[j] += numTemp;
1187 in.close(); m->mothurRemove(tempFile);
1189 for (int i = 0; i < pvalues.size(); i++) { pvalues[i] /= (double)iters; }
1192 pvalues = driver(groupings, groupingsMap, num, indicatorValues, iters);
1193 for (int i = 0; i < pvalues.size(); i++) { pvalues[i] /= (double)iters; }
1198 catch(exception& e) {
1199 m->errorOut(e, "IndicatorCommand", "getPValues");
1204 //**********************************************************************************************************************
1205 //same as above, just data type difference
1206 vector<float> IndicatorCommand::driver(vector< vector<SharedRAbundVector*> >& groupings, map< vector<int>, vector<int> > groupingsMap, int num, vector<float> indicatorValues, int numIters){
1208 vector<float> pvalues;
1209 pvalues.resize(indicatorValues.size(), 0);
1211 for(int i=0;i<numIters;i++){
1212 if (m->control_pressed) { break; }
1213 groupingsMap = randomizeGroupings(groupings, num);
1214 vector<float> randomIndicatorValues = getValues(groupings, groupingsMap);
1216 for (int j = 0; j < indicatorValues.size(); j++) {
1217 if (randomIndicatorValues[j] >= indicatorValues[j]) { pvalues[j]++; }
1223 }catch(exception& e) {
1224 m->errorOut(e, "IndicatorCommand", "driver");
1228 //**********************************************************************************************************************
1229 //same as above, just data type difference
1230 vector<float> IndicatorCommand::getPValues(vector< vector<SharedRAbundVector*> >& groupings, map< vector<int>, vector<int> > groupingsMap, int num, vector<float> indicatorValues){
1232 vector<float> pvalues;
1234 #if defined (__APPLE__) || (__MACH__) || (linux) || (__linux) || (__linux__) || (__unix__) || (__unix)
1235 if(processors == 1){
1236 pvalues = driver(groupings, groupingsMap, num, indicatorValues, iters);
1237 for (int i = 0; i < pvalues.size(); i++) { pvalues[i] /= (double)iters; }
1240 //divide iters between processors
1241 int numItersPerProcessor = iters / processors;
1243 vector<int> processIDS;
1246 //loop through and create all the processes you want
1247 while (process != processors) {
1251 processIDS.push_back(pid); //create map from line number to pid so you can append files in correct order later
1253 }else if (pid == 0){
1254 pvalues = driver(groupings, groupingsMap, num, indicatorValues, numItersPerProcessor);
1256 //pass pvalues to parent
1258 string tempFile = toString(getpid()) + ".pvalues.temp";
1259 m->openOutputFile(tempFile, out);
1262 for (int i = 0; i < pvalues.size(); i++) {
1263 out << pvalues[i] << '\t';
1271 m->mothurOut("[ERROR]: unable to spawn the necessary processes."); m->mothurOutEndLine();
1272 for (int i = 0; i < processIDS.size(); i++) { kill (processIDS[i], SIGINT); }
1278 //special case for last processor in case it doesn't divide evenly
1279 numItersPerProcessor = iters - ((processors-1) * numItersPerProcessor);
1280 pvalues = driver(groupings, groupingsMap, num, indicatorValues, numItersPerProcessor);
1282 //force parent to wait until all the processes are done
1283 for (int i=0;i<processIDS.size();i++) {
1284 int temp = processIDS[i];
1289 for (int i = 0; i < processIDS.size(); i++) {
1291 string tempFile = toString(processIDS[i]) + ".pvalues.temp";
1292 m->openInputFile(tempFile, in);
1294 ////// to do ///////////
1295 int numTemp; numTemp = 0;
1296 for (int j = 0; j < pvalues.size(); j++) {
1297 in >> numTemp; m->gobble(in);
1298 pvalues[j] += numTemp;
1300 in.close(); m->mothurRemove(tempFile);
1302 for (int i = 0; i < pvalues.size(); i++) { pvalues[i] /= (double)iters; }
1305 pvalues = driver(groupings, groupingsMap, num, indicatorValues, iters);
1306 for (int i = 0; i < pvalues.size(); i++) { pvalues[i] /= (double)iters; }
1311 catch(exception& e) {
1312 m->errorOut(e, "IndicatorCommand", "getPValues");
1316 //**********************************************************************************************************************
1317 //swap groups between groupings, in essence randomizing the second column of the design file
1318 map< vector<int>, vector<int> > IndicatorCommand::randomizeGroupings(vector< vector<SharedRAbundVector*> >& groupings, int numLookupGroups){
1321 map< vector<int>, vector<int> > randomGroupings;
1323 for (int i = 0; i < numLookupGroups; i++) {
1324 if (m->control_pressed) {break;}
1326 //get random groups to swap to switch with
1327 //generate random int between 0 and groupings.size()-1
1328 int z = m->getRandomIndex(groupings.size()-1);
1329 int x = m->getRandomIndex(groupings.size()-1);
1330 int a = m->getRandomIndex(groupings[z].size()-1);
1331 int b = m->getRandomIndex(groupings[x].size()-1);
1332 //cout << i << '\t' << z << '\t' << x << '\t' << a << '\t' << b << endl;
1333 //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; }
1338 from.push_back(z); from.push_back(a);
1339 to.push_back(x); to.push_back(b);
1341 randomGroupings[from] = to;
1343 //cout << "done" << endl;
1344 return randomGroupings;
1346 catch(exception& e) {
1347 m->errorOut(e, "IndicatorCommand", "randomizeGroupings");
1351 //**********************************************************************************************************************
1352 //swap groups between groupings, in essence randomizing the second column of the design file
1353 map< vector<int>, vector<int> > IndicatorCommand::randomizeGroupings(vector< vector<SharedRAbundFloatVector*> >& groupings, int numLookupGroups){
1356 map< vector<int>, vector<int> > randomGroupings;
1358 for (int i = 0; i < numLookupGroups; i++) {
1360 //get random groups to swap to switch with
1361 //generate random int between 0 and groupings.size()-1
1362 int z = m->getRandomIndex(groupings.size()-1);
1363 int x = m->getRandomIndex(groupings.size()-1);
1364 int a = m->getRandomIndex(groupings[z].size()-1);
1365 int b = m->getRandomIndex(groupings[x].size()-1);
1366 //cout << i << '\t' << z << '\t' << x << '\t' << a << '\t' << b << endl;
1371 from.push_back(z); from.push_back(a);
1372 to.push_back(x); to.push_back(b);
1374 randomGroupings[from] = to;
1377 return randomGroupings;
1379 catch(exception& e) {
1380 m->errorOut(e, "IndicatorCommand", "randomizeGroupings");
1385 /*****************************************************************/