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 reads a shared or relabund file and a tree or design file, and outputs a .indicator.tre and .indicator.summary file. \n";
42 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";
43 helpString += "The summary file lists the indicator value for each OTU for each node.\n";
44 helpString += "The indicator command parameters are tree, groups, shared, relabund, design and label. The tree parameter is required as well as either shared or relabund.\n";
45 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";
46 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";
47 helpString += "The label parameter indicates at what distance your tree relates to the shared or relabund.\n";
48 helpString += "The processors parameter allows you to specify how many processors you would like to use. The default is 1. \n";
49 helpString += "The iters parameter allows you to set number of randomization for the P value. The default is 1000.";
50 helpString += "The indicator command should be used in the following format: indicator(tree=test.tre, shared=test.shared, label=0.03)\n";
51 helpString += "Note: No spaces between parameter labels (i.e. tree), '=' and parameters (i.e.yourTreefile).\n";
55 m->errorOut(e, "IndicatorCommand", "getHelpString");
60 //**********************************************************************************************************************
61 IndicatorCommand::IndicatorCommand(){
63 abort = true; calledHelp = true;
65 vector<string> tempOutNames;
66 outputTypes["tree"] = tempOutNames;
67 outputTypes["summary"] = tempOutNames;
70 m->errorOut(e, "IndicatorCommand", "IndicatorCommand");
74 //**********************************************************************************************************************
75 IndicatorCommand::IndicatorCommand(string option) {
77 abort = false; calledHelp = false;
79 //allow user to run help
80 if(option == "help") { help(); abort = true; calledHelp = true; }
81 else if(option == "citation") { citation(); abort = true; calledHelp = true;}
84 vector<string> myArray = setParameters();
86 OptionParser parser(option);
87 map<string, string> parameters = parser.getParameters();
89 ValidParameters validParameter;
90 map<string, string>::iterator it;
92 //check to make sure all parameters are valid for command
93 for (it = parameters.begin(); it != parameters.end(); it++) {
94 if (validParameter.isValidParameter(it->first, myArray, it->second) != true) { abort = true; }
99 m->namesOfGroups.clear();
100 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); }
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 convert(temp, iters);
179 temp = validParameter.validFile(parameters, "processors", false); if (temp == "not found"){ temp = m->getProcessors(); }
180 m->setProcessors(temp);
181 convert(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 globaldatas 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
238 SharedUtil* util = new SharedUtil();
239 util->setGroups(Groups, designMap->namesOfGroups);
242 //loop through the Groups and fill Globaldata's Groups with the design file info
243 m->Groups = designMap->getNamesSeqs(Groups);
246 /***************************************************/
247 // use smart distancing to get right sharedRabund //
248 /***************************************************/
249 if (sharedfile != "") {
251 if (m->control_pressed) { if (designfile != "") { delete designMap; } for (int i = 0; i < lookup.size(); i++) { delete lookup[i]; } return 0; }
252 if (lookup[0] == NULL) { m->mothurOut("[ERROR] reading shared file."); m->mothurOutEndLine(); return 0; }
255 if (m->control_pressed) { if (designfile != "") { delete designMap; } for (int i = 0; i < lookupFloat.size(); i++) { delete lookupFloat[i]; } return 0; }
256 if (lookupFloat[0] == NULL) { m->mothurOut("[ERROR] reading relabund file."); m->mothurOutEndLine(); return 0; }
259 //reset Globaldatas groups if needed
260 if (designfile != "") { m->Groups = Groups; }
262 /***************************************************/
263 // reading tree info //
264 /***************************************************/
265 if (treefile != "") {
266 string groupfile = "";
267 m->setTreeFile(treefile);
268 Tree* tree = new Tree(treefile); delete tree; //extracts names from tree to make faked out groupmap
269 treeMap = new TreeMap();
270 bool mismatch = false;
272 for (int i = 0; i < m->Treenames.size(); i++) {
273 //sanity check - is this a group that is not in the sharedfile?
274 if (designfile == "") {
275 if (!(m->inUsersGroups(m->Treenames[i], m->namesOfGroups))) {
276 m->mothurOut("[ERROR]: " + m->Treenames[i] + " is not a group in your shared or relabund file."); m->mothurOutEndLine();
279 treeMap->addSeq(m->Treenames[i], "Group1");
281 vector<string> myGroups; myGroups.push_back(m->Treenames[i]);
282 vector<string> myNames = designMap->getNamesSeqs(myGroups);
284 for(int k = 0; k < myNames.size(); k++) {
285 if (!(m->inUsersGroups(myNames[k], m->namesOfGroups))) {
286 m->mothurOut("[ERROR]: " + myNames[k] + " is not a group in your shared or relabund file."); m->mothurOutEndLine();
290 treeMap->addSeq(m->Treenames[i], "Group1");
294 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; }
296 if (mismatch) { //cleanup and exit
297 if (designfile != "") { delete designMap; }
298 if (sharedfile != "") { for (int i = 0; i < lookup.size(); i++) { delete lookup[i]; } }
299 else { for (int i = 0; i < lookupFloat.size(); i++) { delete lookupFloat[i]; } }
304 read = new ReadNewickTree(treefile);
305 int readOk = read->read(treeMap);
307 if (readOk != 0) { m->mothurOut("Read Terminated."); m->mothurOutEndLine(); delete treeMap; delete read; return 0; }
309 vector<Tree*> T = read->getTrees();
313 if (m->control_pressed) {
314 if (designfile != "") { delete designMap; }
315 if (sharedfile != "") { for (int i = 0; i < lookup.size(); i++) { delete lookup[i]; } }
316 else { for (int i = 0; i < lookupFloat.size(); i++) { delete lookupFloat[i]; } }
317 for (int i = 0; i < T.size(); i++) { delete T[i]; } delete treeMap; return 0;
320 T[0]->assembleTree();
322 /***************************************************/
323 // create ouptut tree - respecting pickedGroups //
324 /***************************************************/
325 Tree* outputTree = new Tree(m->Groups.size(), treeMap);
327 outputTree->getSubTree(T[0], m->Groups);
328 outputTree->assembleTree();
330 //no longer need original tree, we have output tree to use and label
331 for (int i = 0; i < T.size(); i++) { delete T[i]; }
334 if (m->control_pressed) {
335 if (designfile != "") { delete designMap; }
336 if (sharedfile != "") { for (int i = 0; i < lookup.size(); i++) { delete lookup[i]; } }
337 else { for (int i = 0; i < lookupFloat.size(); i++) { delete lookupFloat[i]; } }
338 delete outputTree; delete treeMap; return 0;
341 /***************************************************/
342 // get indicator species values //
343 /***************************************************/
344 GetIndicatorSpecies(outputTree);
345 delete outputTree; delete treeMap;
347 }else { //run with design file only
348 //get indicator species
349 GetIndicatorSpecies();
352 if (designfile != "") { delete designMap; }
353 if (sharedfile != "") { for (int i = 0; i < lookup.size(); i++) { delete lookup[i]; } }
354 else { for (int i = 0; i < lookupFloat.size(); i++) { delete lookupFloat[i]; } }
356 if (m->control_pressed) { for (int i = 0; i < outputNames.size(); i++) { remove(outputNames[i].c_str()); } return 0; }
358 //set tree file as new current treefile
359 if (treefile != "") {
361 itTypes = outputTypes.find("tree");
362 if (itTypes != outputTypes.end()) {
363 if ((itTypes->second).size() != 0) { current = (itTypes->second)[0]; m->setTreeFile(current); }
367 m->mothurOutEndLine(); m->mothurOutEndLine(); m->mothurOut("It took " + toString(time(NULL) - start) + " secs to find the indicator species."); m->mothurOutEndLine();
368 m->mothurOutEndLine();
369 m->mothurOut("Output File Names: "); m->mothurOutEndLine();
370 for (int i = 0; i < outputNames.size(); i++) { m->mothurOut(outputNames[i]); m->mothurOutEndLine(); }
371 m->mothurOutEndLine();
375 catch(exception& e) {
376 m->errorOut(e, "IndicatorCommand", "execute");
380 //**********************************************************************************************************************
381 //divide shared or relabund file by groupings in the design file
382 //report all otu values to file
383 int IndicatorCommand::GetIndicatorSpecies(){
385 string thisOutputDir = outputDir;
386 if (outputDir == "") { thisOutputDir += m->hasPath(inputFileName); }
387 string outputFileName = thisOutputDir + m->getRootName(m->getSimpleName(inputFileName)) + "indicator.summary";
388 outputNames.push_back(outputFileName); outputTypes["summary"].push_back(outputFileName);
391 m->openOutputFile(outputFileName, out);
392 out.setf(ios::fixed, ios::floatfield); out.setf(ios::showpoint);
393 m->mothurOutEndLine(); m->mothurOut("Species\tIndicatorValue\tpValue\n");
396 if (sharedfile != "") { numBins = lookup[0]->getNumBins(); }
397 else { numBins = lookupFloat[0]->getNumBins(); }
399 if (m->control_pressed) { out.close(); return 0; }
401 /*****************************************************/
402 //create vectors containing rabund info //
403 /*****************************************************/
405 vector<float> indicatorValues; //size of numBins
406 vector<float> pValues;
407 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.
409 if (sharedfile != "") {
410 vector< vector<SharedRAbundVector*> > groupings;
411 set<string> groupsAlreadyAdded;
412 vector<SharedRAbundVector*> subset;
415 for (int i = 0; i < designMap->namesOfGroups.size(); i++) {
417 for (int k = 0; k < lookup.size(); k++) {
418 //are you from this grouping?
419 if (designMap->getGroup(lookup[k]->getGroup()) == designMap->namesOfGroups[i]) {
420 subset.push_back(lookup[k]);
421 groupsAlreadyAdded.insert(lookup[k]->getGroup());
424 if (subset.size() != 0) { groupings.push_back(subset); }
428 if (groupsAlreadyAdded.size() != lookup.size()) { m->mothurOut("[ERROR]: could not make proper groupings."); m->mothurOutEndLine(); }
430 indicatorValues = getValues(groupings, randomGroupingsMap);
432 pValues = getPValues(groupings, randomGroupingsMap, lookup.size(), indicatorValues);
434 vector< vector<SharedRAbundFloatVector*> > groupings;
435 set<string> groupsAlreadyAdded;
436 vector<SharedRAbundFloatVector*> subset;
439 for (int i = 0; i < designMap->namesOfGroups.size(); i++) {
440 for (int k = 0; k < lookupFloat.size(); k++) {
441 //are you from this grouping?
442 if (designMap->getGroup(lookupFloat[k]->getGroup()) == designMap->namesOfGroups[i]) {
443 subset.push_back(lookupFloat[k]);
444 groupsAlreadyAdded.insert(lookupFloat[k]->getGroup());
447 if (subset.size() != 0) { groupings.push_back(subset); }
451 if (groupsAlreadyAdded.size() != lookupFloat.size()) { m->mothurOut("[ERROR]: could not make proper groupings."); m->mothurOutEndLine(); }
453 indicatorValues = getValues(groupings, randomGroupingsMap);
455 pValues = getPValues(groupings, randomGroupingsMap, lookupFloat.size(), indicatorValues);
458 if (m->control_pressed) { out.close(); return 0; }
461 /******************************************************/
462 //output indicator values to table form //
463 /*****************************************************/
464 out << "OTU\tIndicator_Value\tpValue" << endl;
465 for (int j = 0; j < indicatorValues.size(); j++) {
467 if (m->control_pressed) { out.close(); return 0; }
469 out << (j+1) << '\t' << indicatorValues[j] << '\t';
471 if (pValues[j] > (1/(float)iters)) { out << pValues[j] << endl; }
472 else { out << "<" << (1/(float)iters) << endl; }
474 if (pValues[j] <= 0.05) {
475 cout << "OTU" << j+1 << '\t' << indicatorValues[j] << '\t';
476 string pValueString = "<" + toString((1/(float)iters));
477 if (pValues[j] > (1/(float)iters)) { pValueString = toString(pValues[j]); cout << pValues[j];}
478 else { cout << "<" << (1/(float)iters); }
479 m->mothurOutJustToLog("OTU" + toString(j+1) + "\t" + toString(indicatorValues[j]) + "\t" + pValueString);
480 m->mothurOutEndLine();
488 catch(exception& e) {
489 m->errorOut(e, "IndicatorCommand", "GetIndicatorSpecies");
493 //**********************************************************************************************************************
494 //traverse tree finding indicator species values for each otu at each node
495 //label node with otu number that has highest indicator value
496 //report all otu values to file
497 int IndicatorCommand::GetIndicatorSpecies(Tree*& T){
500 string thisOutputDir = outputDir;
501 if (outputDir == "") { thisOutputDir += m->hasPath(inputFileName); }
502 string outputFileName = thisOutputDir + m->getRootName(m->getSimpleName(inputFileName)) + "indicator.summary";
503 outputNames.push_back(outputFileName); outputTypes["summary"].push_back(outputFileName);
506 m->openOutputFile(outputFileName, out);
507 out.setf(ios::fixed, ios::floatfield); out.setf(ios::showpoint);
510 if (sharedfile != "") { numBins = lookup[0]->getNumBins(); }
511 else { numBins = lookupFloat[0]->getNumBins(); }
515 for (int i = 0; i < numBins; i++) { out << "OTU" << (i+1) << "_IndValue" << '\t' << "pValue" << '\t'; }
518 m->mothurOutEndLine(); m->mothurOut("Node\tSpecies\tIndicatorValue\tpValue\n");
520 string treeOutputDir = outputDir;
521 if (outputDir == "") { treeOutputDir += m->hasPath(treefile); }
522 string outputTreeFileName = treeOutputDir + m->getRootName(m->getSimpleName(treefile)) + "indicator.tre";
525 //create a map from tree node index to names of descendants, save time later to know which sharedRabund you need
526 map<int, set<string> > nodeToDescendants;
527 map<int, set<int> > descendantNodes;
528 for (int i = 0; i < T->getNumNodes(); i++) {
529 if (m->control_pressed) { return 0; }
531 nodeToDescendants[i] = getDescendantList(T, i, nodeToDescendants, descendantNodes);
534 //you need the distances to leaf to decide grouping below
535 //this will also set branch lengths if the tree does not include them
536 map<int, float> distToRoot = getDistToRoot(T);
539 for (int i = T->getNumLeaves(); i < T->getNumNodes(); i++) {
540 //cout << endl << i+1 << endl;
541 if (m->control_pressed) { out.close(); return 0; }
543 /*****************************************************/
544 //create vectors containing rabund info //
545 /*****************************************************/
547 vector<float> indicatorValues; //size of numBins
548 vector<float> pValues;
549 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.
551 if (sharedfile != "") {
552 vector< vector<SharedRAbundVector*> > groupings;
554 //get nodes that will be a valid grouping
555 //you are valid if you are not one of my descendants
556 //AND your distToRoot is >= mine
557 //AND you were not added as part of a larger grouping. Largest nodes are added first.
559 set<string> groupsAlreadyAdded;
560 //create a grouping with my grouping
561 vector<SharedRAbundVector*> subset;
563 int doneCount = nodeToDescendants[i].size();
564 for (int k = 0; k < lookup.size(); k++) {
565 //is this descendant of i
566 if ((nodeToDescendants[i].count(lookup[k]->getGroup()) != 0)) {
567 subset.push_back(lookup[k]);
568 groupsAlreadyAdded.insert(lookup[k]->getGroup());
571 if (count == doneCount) { break; } //quit once you get the rabunds for this grouping
573 if (subset.size() != 0) { groupings.push_back(subset); }
576 for (int j = (T->getNumNodes()-1); j >= 0; j--) {
579 if ((descendantNodes[i].count(j) == 0) && (distToRoot[j] >= distToRoot[i])) {
580 vector<SharedRAbundVector*> subset;
582 int doneCount = nodeToDescendants[j].size();
583 for (int k = 0; k < lookup.size(); k++) {
584 //is this descendant of j, and we didn't already add this as part of a larger grouping
585 if ((nodeToDescendants[j].count(lookup[k]->getGroup()) != 0) && (groupsAlreadyAdded.count(lookup[k]->getGroup()) == 0)) {
586 subset.push_back(lookup[k]);
587 groupsAlreadyAdded.insert(lookup[k]->getGroup());
590 if (count == doneCount) { break; } //quit once you get the rabunds for this grouping
593 //if subset.size == 0 then the node was added as part of a larger grouping
594 if (subset.size() != 0) { groupings.push_back(subset); }
598 if (groupsAlreadyAdded.size() != lookup.size()) { m->mothurOut("[ERROR]: could not make proper groupings."); m->mothurOutEndLine(); }
600 indicatorValues = getValues(groupings, randomGroupingsMap);
602 pValues = getPValues(groupings, randomGroupingsMap, lookup.size(), indicatorValues);
604 vector< vector<SharedRAbundFloatVector*> > groupings;
606 //get nodes that will be a valid grouping
607 //you are valid if you are not one of my descendants
608 //AND your distToRoot is >= mine
609 //AND you were not added as part of a larger grouping. Largest nodes are added first.
611 set<string> groupsAlreadyAdded;
612 //create a grouping with my grouping
613 vector<SharedRAbundFloatVector*> subset;
615 int doneCount = nodeToDescendants[i].size();
616 for (int k = 0; k < lookupFloat.size(); k++) {
617 //is this descendant of i
618 if ((nodeToDescendants[i].count(lookupFloat[k]->getGroup()) != 0)) {
619 subset.push_back(lookupFloat[k]);
620 groupsAlreadyAdded.insert(lookupFloat[k]->getGroup());
623 if (count == doneCount) { break; } //quit once you get the rabunds for this grouping
625 if (subset.size() != 0) { groupings.push_back(subset); }
627 for (int j = (T->getNumNodes()-1); j >= 0; j--) {
628 if ((descendantNodes[i].count(j) == 0) && (distToRoot[j] >= distToRoot[i])) {
629 vector<SharedRAbundFloatVector*> subset;
631 int doneCount = nodeToDescendants[j].size();
632 for (int k = 0; k < lookupFloat.size(); k++) {
633 //is this descendant of j, and we didn't already add this as part of a larger grouping
634 if ((nodeToDescendants[j].count(lookupFloat[k]->getGroup()) != 0) && (groupsAlreadyAdded.count(lookupFloat[k]->getGroup()) == 0)) {
635 subset.push_back(lookupFloat[k]);
636 groupsAlreadyAdded.insert(lookupFloat[k]->getGroup());
639 if (count == doneCount) { break; } //quit once you get the rabunds for this grouping
642 //if subset.size == 0 then the node was added as part of a larger grouping
643 if (subset.size() != 0) { groupings.push_back(subset); }
647 if (groupsAlreadyAdded.size() != lookupFloat.size()) { m->mothurOut("[ERROR]: could not make proper groupings."); m->mothurOutEndLine(); }
649 indicatorValues = getValues(groupings, randomGroupingsMap);
651 pValues = getPValues(groupings, randomGroupingsMap, lookupFloat.size(), indicatorValues);
654 if (m->control_pressed) { out.close(); return 0; }
657 /******************************************************/
658 //output indicator values to table form + label tree //
659 /*****************************************************/
660 out << (i+1) << '\t';
661 for (int j = 0; j < indicatorValues.size(); j++) {
663 if (m->control_pressed) { out.close(); return 0; }
665 if (pValues[j] < (1/(float)iters)) {
666 out << indicatorValues[j] << '\t' << '<' << (1/(float)iters) << '\t';
668 out << indicatorValues[j] << '\t' << pValues[j] << '\t';
671 if (pValues[j] <= 0.05) {
672 cout << i+1 << "\tOTU" << j+1 << '\t' << indicatorValues[j] << '\t';
673 string pValueString = "<" + toString((1/(float)iters));
674 if (pValues[j] > (1/(float)iters)) { pValueString = toString(pValues[j]); cout << pValues[j];}
675 else { cout << "<" << (1/(float)iters); }
676 m->mothurOutJustToLog(toString(i) + "\tOTU" + toString(j+1) + "\t" + toString(indicatorValues[j]) + "\t" + pValueString);
677 m->mothurOutEndLine();
682 T->tree[i].setLabel((i+1));
687 m->openOutputFile(outputTreeFileName, outTree);
688 outputNames.push_back(outputTreeFileName); outputTypes["tree"].push_back(outputTreeFileName);
690 T->print(outTree, "both");
695 catch(exception& e) {
696 m->errorOut(e, "IndicatorCommand", "GetIndicatorSpecies");
700 //**********************************************************************************************************************
701 vector<float> IndicatorCommand::getValues(vector< vector<SharedRAbundFloatVector*> >& groupings, map< vector<int>, vector<int> > groupingsMap){
703 vector<float> values;
704 map< vector<int>, vector<int> >::iterator it;
707 for (int i = 0; i < groupings[0][0]->getNumBins(); i++) {
709 if (m->control_pressed) { return values; }
712 float AijDenominator = 0.0;
715 //get overall abundance of each grouping
716 for (int j = 0; j < groupings.size(); j++) {
718 float totalAbund = 0;
720 for (int k = 0; k < groupings[j].size(); k++) {
721 vector<int> temp; temp.push_back(j); temp.push_back(k);
722 it = groupingsMap.find(temp);
724 if (it == groupingsMap.end()) { //this one didnt get moved
725 totalAbund += groupings[j][k]->getAbundance(i);
726 if (groupings[j][k]->getAbundance(i) != 0.0) { numNotZero++; }
728 totalAbund += groupings[(it->second)[0]][(it->second)[1]]->getAbundance(i);
729 if (groupings[(it->second)[0]][(it->second)[1]]->getAbundance(i) != 0.0) { numNotZero++; }
735 float Aij = (totalAbund / (float) groupings[j].size());
736 terms.push_back(Aij);
738 //percentage of sites represented
739 Bij.push_back(numNotZero / (float) groupings[j].size());
741 AijDenominator += Aij;
744 float maxIndVal = 0.0;
745 for (int j = 0; j < terms.size(); j++) {
746 float thisAij = (terms[j] / AijDenominator); //relative abundance
747 float thisValue = thisAij * Bij[j] * 100.0;
750 if (thisValue > maxIndVal) { maxIndVal = thisValue; }
753 values.push_back(maxIndVal);
758 catch(exception& e) {
759 m->errorOut(e, "IndicatorCommand", "getValues");
763 //**********************************************************************************************************************
764 //same as above, just data type difference
765 vector<float> IndicatorCommand::getValues(vector< vector<SharedRAbundVector*> >& groupings, map< vector<int>, vector<int> > groupingsMap){
767 vector<float> values;
769 /*for (int j = 0; j < groupings.size(); j++) {
770 cout << "grouping " << j << endl;
771 for (int k = 0; k < groupings[j].size(); k++) {
772 cout << groupings[j][k]->getGroup() << endl;
775 map< vector<int>, vector<int> >::iterator it;
778 for (int i = 0; i < groupings[0][0]->getNumBins(); i++) {
780 float AijDenominator = 0.0;
782 //get overall abundance of each grouping
783 for (int j = 0; j < groupings.size(); j++) {
785 int totalAbund = 0.0;
787 for (int k = 0; k < groupings[j].size(); k++) {
788 vector<int> temp; temp.push_back(j); temp.push_back(k);
789 it = groupingsMap.find(temp);
791 if (it == groupingsMap.end()) { //this one didnt get moved
792 totalAbund += groupings[j][k]->getAbundance(i);
793 if (groupings[j][k]->getAbundance(i) != 0.0) { numNotZero++; }
795 totalAbund += groupings[(it->second)[0]][(it->second)[1]]->getAbundance(i);
796 if (groupings[(it->second)[0]][(it->second)[1]]->getAbundance(i) != 0.0) { numNotZero++; }
802 float Aij = (totalAbund / (float) groupings[j].size());
803 terms.push_back(Aij);
805 //percentage of sites represented
806 Bij.push_back(numNotZero / (float) groupings[j].size());
808 AijDenominator += Aij;
811 float maxIndVal = 0.0;
812 for (int j = 0; j < terms.size(); j++) {
813 float thisAij = (terms[j] / AijDenominator); //relative abundance
814 float thisValue = thisAij * Bij[j] * 100.0;
817 if (thisValue > maxIndVal) { maxIndVal = thisValue; }
820 values.push_back(maxIndVal);
825 catch(exception& e) {
826 m->errorOut(e, "IndicatorCommand", "getValues");
830 //**********************************************************************************************************************
831 //you need the distances to root to decide groupings
832 //this will also set branch lengths if the tree does not include them
833 map<int, float> IndicatorCommand::getDistToRoot(Tree*& T){
835 map<int, float> dists;
837 bool hasBranchLengths = false;
838 for (int i = 0; i < T->getNumNodes(); i++) {
839 if (T->tree[i].getBranchLength() > 0.0) { hasBranchLengths = true; break; }
842 //set branchlengths if needed
843 if (!hasBranchLengths) {
844 for (int i = 0; i < T->getNumNodes(); i++) {
846 int lc = T->tree[i].getLChild();
847 int rc = T->tree[i].getRChild();
849 if (lc == -1) { // you are a leaf
850 //if you are a leaf set you priliminary length to 1.0, this may adjust later
851 T->tree[i].setBranchLength(1.0);
853 }else{ // you are an internal node
854 //look at your children's length to leaf
855 float ldist = dists[lc];
856 float rdist = dists[rc];
858 float greater = ldist;
859 if (rdist > greater) { greater = rdist; dists[i] = ldist + 1.0;}
860 else { dists[i] = rdist + 1.0; }
863 //branch length = difference + 1
864 T->tree[lc].setBranchLength((abs(ldist-greater) + 1.0));
865 T->tree[rc].setBranchLength((abs(rdist-greater) + 1.0));
872 for (int i = 0; i < T->getNumNodes(); i++) {
877 while(T->tree[index].getParent() != -1){
878 if (T->tree[index].getBranchLength() != -1) {
879 sum += abs(T->tree[index].getBranchLength());
881 index = T->tree[index].getParent();
889 catch(exception& e) {
890 m->errorOut(e, "IndicatorCommand", "getLengthToLeaf");
894 //**********************************************************************************************************************
895 set<string> IndicatorCommand::getDescendantList(Tree*& T, int i, map<int, set<string> > descendants, map<int, set<int> >& nodes){
899 set<string>::iterator it;
901 int lc = T->tree[i].getLChild();
902 int rc = T->tree[i].getRChild();
904 if (lc == -1) { //you are a leaf your only descendant is yourself
905 set<int> temp; temp.insert(i);
908 if (designfile == "") {
909 names.insert(T->tree[i].getName());
911 vector<string> myGroup; myGroup.push_back(T->tree[i].getName());
912 vector<string> myReps = designMap->getNamesSeqs(myGroup);
913 for (int k = 0; k < myReps.size(); k++) {
914 names.insert(myReps[k]);
918 }else{ //your descedants are the combination of your childrens descendants
919 names = descendants[lc];
920 nodes[i] = nodes[lc];
921 for (it = descendants[rc].begin(); it != descendants[rc].end(); it++) {
924 for (set<int>::iterator itNum = nodes[rc].begin(); itNum != nodes[rc].end(); itNum++) {
925 nodes[i].insert(*itNum);
927 //you are your own descendant
933 catch(exception& e) {
934 m->errorOut(e, "IndicatorCommand", "getDescendantList");
938 //**********************************************************************************************************************
939 int IndicatorCommand::getShared(){
941 InputData* input = new InputData(sharedfile, "sharedfile");
942 lookup = input->getSharedRAbundVectors();
943 string lastLabel = lookup[0]->getLabel();
945 if (label == "") { label = lastLabel; delete input; return 0; }
947 //if the users enters label "0.06" and there is no "0.06" in their file use the next lowest label.
948 set<string> labels; labels.insert(label);
949 set<string> processedLabels;
950 set<string> userLabels = labels;
952 //as long as you are not at the end of the file or done wih the lines you want
953 while((lookup[0] != NULL) && (userLabels.size() != 0)) {
954 if (m->control_pressed) { delete input; return 0; }
956 if(labels.count(lookup[0]->getLabel()) == 1){
957 processedLabels.insert(lookup[0]->getLabel());
958 userLabels.erase(lookup[0]->getLabel());
962 if ((m->anyLabelsToProcess(lookup[0]->getLabel(), userLabels, "") == true) && (processedLabels.count(lastLabel) != 1)) {
963 string saveLabel = lookup[0]->getLabel();
965 for (int i = 0; i < lookup.size(); i++) { delete lookup[i]; }
966 lookup = input->getSharedRAbundVectors(lastLabel);
968 processedLabels.insert(lookup[0]->getLabel());
969 userLabels.erase(lookup[0]->getLabel());
971 //restore real lastlabel to save below
972 lookup[0]->setLabel(saveLabel);
976 lastLabel = lookup[0]->getLabel();
978 //get next line to process
979 //prevent memory leak
980 for (int i = 0; i < lookup.size(); i++) { delete lookup[i]; }
981 lookup = input->getSharedRAbundVectors();
985 if (m->control_pressed) { delete input; return 0; }
987 //output error messages about any remaining user labels
988 set<string>::iterator it;
989 bool needToRun = false;
990 for (it = userLabels.begin(); it != userLabels.end(); it++) {
991 m->mothurOut("Your file does not include the label " + *it);
992 if (processedLabels.count(lastLabel) != 1) {
993 m->mothurOut(". I will use " + lastLabel + "."); m->mothurOutEndLine();
996 m->mothurOut(". Please refer to " + lastLabel + "."); m->mothurOutEndLine();
1000 //run last label if you need to
1001 if (needToRun == true) {
1002 for (int i = 0; i < lookup.size(); i++) { if (lookup[i] != NULL) { delete lookup[i]; } }
1003 lookup = input->getSharedRAbundVectors(lastLabel);
1009 catch(exception& e) {
1010 m->errorOut(e, "IndicatorCommand", "getShared");
1014 //**********************************************************************************************************************
1015 int IndicatorCommand::getSharedFloat(){
1017 InputData* input = new InputData(relabundfile, "relabund");
1018 lookupFloat = input->getSharedRAbundFloatVectors();
1019 string lastLabel = lookupFloat[0]->getLabel();
1021 if (label == "") { label = lastLabel; delete input; return 0; }
1023 //if the users enters label "0.06" and there is no "0.06" in their file use the next lowest label.
1024 set<string> labels; labels.insert(label);
1025 set<string> processedLabels;
1026 set<string> userLabels = labels;
1028 //as long as you are not at the end of the file or done wih the lines you want
1029 while((lookupFloat[0] != NULL) && (userLabels.size() != 0)) {
1031 if (m->control_pressed) { delete input; return 0; }
1033 if(labels.count(lookupFloat[0]->getLabel()) == 1){
1034 processedLabels.insert(lookupFloat[0]->getLabel());
1035 userLabels.erase(lookupFloat[0]->getLabel());
1039 if ((m->anyLabelsToProcess(lookupFloat[0]->getLabel(), userLabels, "") == true) && (processedLabels.count(lastLabel) != 1)) {
1040 string saveLabel = lookupFloat[0]->getLabel();
1042 for (int i = 0; i < lookupFloat.size(); i++) { delete lookupFloat[i]; }
1043 lookupFloat = input->getSharedRAbundFloatVectors(lastLabel);
1045 processedLabels.insert(lookupFloat[0]->getLabel());
1046 userLabels.erase(lookupFloat[0]->getLabel());
1048 //restore real lastlabel to save below
1049 lookupFloat[0]->setLabel(saveLabel);
1053 lastLabel = lookupFloat[0]->getLabel();
1055 //get next line to process
1056 //prevent memory leak
1057 for (int i = 0; i < lookupFloat.size(); i++) { delete lookupFloat[i]; }
1058 lookupFloat = input->getSharedRAbundFloatVectors();
1062 if (m->control_pressed) { delete input; return 0; }
1064 //output error messages about any remaining user labels
1065 set<string>::iterator it;
1066 bool needToRun = false;
1067 for (it = userLabels.begin(); it != userLabels.end(); it++) {
1068 m->mothurOut("Your file does not include the label " + *it);
1069 if (processedLabels.count(lastLabel) != 1) {
1070 m->mothurOut(". I will use " + lastLabel + "."); m->mothurOutEndLine();
1073 m->mothurOut(". Please refer to " + lastLabel + "."); m->mothurOutEndLine();
1077 //run last label if you need to
1078 if (needToRun == true) {
1079 for (int i = 0; i < lookupFloat.size(); i++) { if (lookupFloat[i] != NULL) { delete lookupFloat[i]; } }
1080 lookupFloat = input->getSharedRAbundFloatVectors(lastLabel);
1086 catch(exception& e) {
1087 m->errorOut(e, "IndicatorCommand", "getShared");
1091 //**********************************************************************************************************************
1092 vector<float> IndicatorCommand::driver(vector< vector<SharedRAbundFloatVector*> >& groupings, map< vector<int>, vector<int> > groupingsMap, int num, vector<float> indicatorValues, int numIters){
1094 vector<float> pvalues;
1095 pvalues.resize(indicatorValues.size(), 0);
1097 for(int i=0;i<numIters;i++){
1098 if (m->control_pressed) { break; }
1099 groupingsMap = randomizeGroupings(groupings, num);
1100 vector<float> randomIndicatorValues = getValues(groupings, groupingsMap);
1102 for (int j = 0; j < indicatorValues.size(); j++) {
1103 if (randomIndicatorValues[j] >= indicatorValues[j]) { pvalues[j]++; }
1109 }catch(exception& e) {
1110 m->errorOut(e, "IndicatorCommand", "driver");
1114 //**********************************************************************************************************************
1115 vector<float> IndicatorCommand::getPValues(vector< vector<SharedRAbundFloatVector*> >& groupings, map< vector<int>, vector<int> > groupingsMap, int num, vector<float> indicatorValues){
1117 vector<float> pvalues;
1119 #if defined (__APPLE__) || (__MACH__) || (linux) || (__linux)
1120 if(processors == 1){
1121 pvalues = driver(groupings, groupingsMap, num, indicatorValues, iters);
1122 for (int i = 0; i < pvalues.size(); i++) { pvalues[i] /= (double)iters; }
1125 //divide iters between processors
1126 int numItersPerProcessor = iters / processors;
1128 vector<int> processIDS;
1132 //loop through and create all the processes you want
1133 while (process != processors) {
1137 processIDS.push_back(pid); //create map from line number to pid so you can append files in correct order later
1139 }else if (pid == 0){
1140 pvalues = driver(groupings, groupingsMap, num, indicatorValues, numItersPerProcessor);
1142 //pass pvalues to parent
1144 string tempFile = toString(getpid()) + ".pvalues.temp";
1145 m->openOutputFile(tempFile, out);
1148 for (int i = 0; i < pvalues.size(); i++) {
1149 out << pvalues[i] << '\t';
1157 m->mothurOut("[ERROR]: unable to spawn the necessary processes."); m->mothurOutEndLine();
1158 for (int i = 0; i < processIDS.size(); i++) { kill (processIDS[i], SIGINT); }
1164 //special case for last processor in case it doesn't divide evenly
1165 numItersPerProcessor = iters - ((processors-1) * numItersPerProcessor);
1166 pvalues = driver(groupings, groupingsMap, num, indicatorValues, numItersPerProcessor);
1168 //force parent to wait until all the processes are done
1169 for (int i=0;i<processIDS.size();i++) {
1170 int temp = processIDS[i];
1175 for (int i = 0; i < processIDS.size(); i++) {
1177 string tempFile = toString(processIDS[i]) + ".pvalues.temp";
1178 m->openInputFile(tempFile, in);
1180 ////// to do ///////////
1181 int numTemp; numTemp = 0;
1182 for (int j = 0; j < pvalues.size(); j++) {
1183 in >> numTemp; m->gobble(in);
1184 pvalues[j] += numTemp;
1186 in.close(); remove(tempFile.c_str());
1188 for (int i = 0; i < pvalues.size(); i++) { pvalues[i] /= (double)iters; }
1191 pvalues = driver(groupings, groupingsMap, num, indicatorValues, iters);
1192 for (int i = 0; i < pvalues.size(); i++) { pvalues[i] /= (double)iters; }
1197 catch(exception& e) {
1198 m->errorOut(e, "IndicatorCommand", "getPValues");
1203 //**********************************************************************************************************************
1204 //same as above, just data type difference
1205 vector<float> IndicatorCommand::driver(vector< vector<SharedRAbundVector*> >& groupings, map< vector<int>, vector<int> > groupingsMap, int num, vector<float> indicatorValues, int numIters){
1207 vector<float> pvalues;
1208 pvalues.resize(indicatorValues.size(), 0);
1210 for(int i=0;i<numIters;i++){
1211 if (m->control_pressed) { break; }
1212 groupingsMap = randomizeGroupings(groupings, num);
1213 vector<float> randomIndicatorValues = getValues(groupings, groupingsMap);
1215 for (int j = 0; j < indicatorValues.size(); j++) {
1216 if (randomIndicatorValues[j] >= indicatorValues[j]) { pvalues[j]++; }
1222 }catch(exception& e) {
1223 m->errorOut(e, "IndicatorCommand", "driver");
1227 //**********************************************************************************************************************
1228 //same as above, just data type difference
1229 vector<float> IndicatorCommand::getPValues(vector< vector<SharedRAbundVector*> >& groupings, map< vector<int>, vector<int> > groupingsMap, int num, vector<float> indicatorValues){
1231 vector<float> pvalues;
1233 #if defined (__APPLE__) || (__MACH__) || (linux) || (__linux)
1234 if(processors == 1){
1235 pvalues = driver(groupings, groupingsMap, num, indicatorValues, iters);
1236 for (int i = 0; i < pvalues.size(); i++) { pvalues[i] /= (double)iters; }
1239 //divide iters between processors
1240 int numItersPerProcessor = iters / processors;
1242 vector<int> processIDS;
1245 //loop through and create all the processes you want
1246 while (process != processors) {
1250 processIDS.push_back(pid); //create map from line number to pid so you can append files in correct order later
1252 }else if (pid == 0){
1253 pvalues = driver(groupings, groupingsMap, num, indicatorValues, numItersPerProcessor);
1255 //pass pvalues to parent
1257 string tempFile = toString(getpid()) + ".pvalues.temp";
1258 m->openOutputFile(tempFile, out);
1261 for (int i = 0; i < pvalues.size(); i++) {
1262 out << pvalues[i] << '\t';
1270 m->mothurOut("[ERROR]: unable to spawn the necessary processes."); m->mothurOutEndLine();
1271 for (int i = 0; i < processIDS.size(); i++) { kill (processIDS[i], SIGINT); }
1277 //special case for last processor in case it doesn't divide evenly
1278 numItersPerProcessor = iters - ((processors-1) * numItersPerProcessor);
1279 pvalues = driver(groupings, groupingsMap, num, indicatorValues, numItersPerProcessor);
1281 //force parent to wait until all the processes are done
1282 for (int i=0;i<processIDS.size();i++) {
1283 int temp = processIDS[i];
1288 for (int i = 0; i < processIDS.size(); i++) {
1290 string tempFile = toString(processIDS[i]) + ".pvalues.temp";
1291 m->openInputFile(tempFile, in);
1293 ////// to do ///////////
1294 int numTemp; numTemp = 0;
1295 for (int j = 0; j < pvalues.size(); j++) {
1296 in >> numTemp; m->gobble(in);
1297 pvalues[j] += numTemp;
1299 in.close(); remove(tempFile.c_str());
1301 for (int i = 0; i < pvalues.size(); i++) { pvalues[i] /= (double)iters; }
1304 pvalues = driver(groupings, groupingsMap, num, indicatorValues, iters);
1305 for (int i = 0; i < pvalues.size(); i++) { pvalues[i] /= (double)iters; }
1310 catch(exception& e) {
1311 m->errorOut(e, "IndicatorCommand", "getPValues");
1315 //**********************************************************************************************************************
1316 //swap groups between groupings, in essence randomizing the second column of the design file
1317 map< vector<int>, vector<int> > IndicatorCommand::randomizeGroupings(vector< vector<SharedRAbundVector*> >& groupings, int numLookupGroups){
1320 map< vector<int>, vector<int> > randomGroupings;
1322 for (int i = 0; i < numLookupGroups; i++) {
1323 if (m->control_pressed) {break;}
1325 //get random groups to swap to switch with
1326 //generate random int between 0 and groupings.size()-1
1327 int z = m->getRandomIndex(groupings.size()-1);
1328 int x = m->getRandomIndex(groupings.size()-1);
1329 int a = m->getRandomIndex(groupings[z].size()-1);
1330 int b = m->getRandomIndex(groupings[x].size()-1);
1331 //cout << i << '\t' << z << '\t' << x << '\t' << a << '\t' << b << endl;
1332 //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; }
1337 from.push_back(z); from.push_back(a);
1338 to.push_back(x); to.push_back(b);
1340 randomGroupings[from] = to;
1342 //cout << "done" << endl;
1343 return randomGroupings;
1345 catch(exception& e) {
1346 m->errorOut(e, "IndicatorCommand", "randomizeGroupings");
1350 //**********************************************************************************************************************
1351 //swap groups between groupings, in essence randomizing the second column of the design file
1352 map< vector<int>, vector<int> > IndicatorCommand::randomizeGroupings(vector< vector<SharedRAbundFloatVector*> >& groupings, int numLookupGroups){
1355 map< vector<int>, vector<int> > randomGroupings;
1357 for (int i = 0; i < numLookupGroups; i++) {
1359 //get random groups to swap to switch with
1360 //generate random int between 0 and groupings.size()-1
1361 int z = m->getRandomIndex(groupings.size()-1);
1362 int x = m->getRandomIndex(groupings.size()-1);
1363 int a = m->getRandomIndex(groupings[z].size()-1);
1364 int b = m->getRandomIndex(groupings[x].size()-1);
1365 //cout << i << '\t' << z << '\t' << x << '\t' << a << '\t' << b << endl;
1370 from.push_back(z); from.push_back(a);
1371 to.push_back(x); to.push_back(b);
1373 randomGroupings[from] = to;
1376 return randomGroupings;
1378 catch(exception& e) {
1379 m->errorOut(e, "IndicatorCommand", "randomizeGroupings");
1384 /*****************************************************************/