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");
60 //**********************************************************************************************************************
61 string IndicatorCommand::getOutputFileNameTag(string type, string inputName=""){
63 string outputFileName = "";
64 map<string, vector<string> >::iterator it;
66 //is this a type this command creates
67 it = outputTypes.find(type);
68 if (it == outputTypes.end()) { m->mothurOut("[ERROR]: this command doesn't create a " + type + " output file.\n"); }
70 if (type == "tree") { outputFileName = "indicator.tre"; }
71 else if (type == "summary") { outputFileName = "indicator.summary"; }
72 else { m->mothurOut("[ERROR]: No definition for type " + type + " output file tag.\n"); m->control_pressed = true; }
74 return outputFileName;
77 m->errorOut(e, "IndicatorCommand", "getOutputFileNameTag");
81 //**********************************************************************************************************************
82 IndicatorCommand::IndicatorCommand(){
84 abort = true; calledHelp = true;
86 vector<string> tempOutNames;
87 outputTypes["tree"] = tempOutNames;
88 outputTypes["summary"] = tempOutNames;
91 m->errorOut(e, "IndicatorCommand", "IndicatorCommand");
95 //**********************************************************************************************************************
96 IndicatorCommand::IndicatorCommand(string option) {
98 abort = false; calledHelp = false;
100 //allow user to run help
101 if(option == "help") { help(); abort = true; calledHelp = true; }
102 else if(option == "citation") { citation(); abort = true; calledHelp = true;}
105 vector<string> myArray = setParameters();
107 OptionParser parser(option);
108 map<string, string> parameters = parser.getParameters();
110 ValidParameters validParameter;
111 map<string, string>::iterator it;
113 //check to make sure all parameters are valid for command
114 for (it = parameters.begin(); it != parameters.end(); it++) {
115 if (validParameter.isValidParameter(it->first, myArray, it->second) != true) { abort = true; }
121 m->Treenames.clear();
123 vector<string> tempOutNames;
124 outputTypes["tree"] = tempOutNames;
125 outputTypes["summary"] = tempOutNames;
127 //if the user changes the input directory command factory will send this info to us in the output parameter
128 string inputDir = validParameter.validFile(parameters, "inputdir", false);
129 if (inputDir == "not found"){ inputDir = ""; }
132 it = parameters.find("tree");
133 //user has given a template file
134 if(it != parameters.end()){
135 path = m->hasPath(it->second);
136 //if the user has not given a path then, add inputdir. else leave path alone.
137 if (path == "") { parameters["tree"] = inputDir + it->second; }
140 it = parameters.find("shared");
141 //user has given a template file
142 if(it != parameters.end()){
143 path = m->hasPath(it->second);
144 //if the user has not given a path then, add inputdir. else leave path alone.
145 if (path == "") { parameters["shared"] = inputDir + it->second; }
148 it = parameters.find("relabund");
149 //user has given a template file
150 if(it != parameters.end()){
151 path = m->hasPath(it->second);
152 //if the user has not given a path then, add inputdir. else leave path alone.
153 if (path == "") { parameters["relabund"] = inputDir + it->second; }
156 it = parameters.find("design");
157 //user has given a template file
158 if(it != parameters.end()){
159 path = m->hasPath(it->second);
160 //if the user has not given a path then, add inputdir. else leave path alone.
161 if (path == "") { parameters["design"] = inputDir + it->second; }
165 outputDir = validParameter.validFile(parameters, "outputdir", false); if (outputDir == "not found"){ outputDir = ""; }
167 //check for required parameters
168 treefile = validParameter.validFile(parameters, "tree", true);
169 if (treefile == "not open") { treefile = ""; abort = true; }
170 else if (treefile == "not found") { treefile = ""; }
171 else { m->setTreeFile(treefile); }
173 sharedfile = validParameter.validFile(parameters, "shared", true);
174 if (sharedfile == "not open") { abort = true; }
175 else if (sharedfile == "not found") { sharedfile = ""; }
176 else { inputFileName = sharedfile; m->setSharedFile(sharedfile); }
178 relabundfile = validParameter.validFile(parameters, "relabund", true);
179 if (relabundfile == "not open") { abort = true; }
180 else if (relabundfile == "not found") { relabundfile = ""; }
181 else { inputFileName = relabundfile; m->setRelAbundFile(relabundfile); }
183 designfile = validParameter.validFile(parameters, "design", true);
184 if (designfile == "not open") { designfile = ""; abort = true; }
185 else if (designfile == "not found") { designfile = ""; }
186 else { m->setDesignFile(designfile); }
188 groups = validParameter.validFile(parameters, "groups", false);
189 if (groups == "not found") { groups = ""; Groups.push_back("all"); }
190 else { m->splitAtDash(groups, Groups); }
191 m->setGroups(Groups);
193 label = validParameter.validFile(parameters, "label", false);
194 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=""; }
196 string temp = validParameter.validFile(parameters, "iters", false); if (temp == "not found") { temp = "1000"; }
197 m->mothurConvert(temp, iters);
199 temp = validParameter.validFile(parameters, "processors", false); if (temp == "not found"){ temp = m->getProcessors(); }
200 m->setProcessors(temp);
201 m->mothurConvert(temp, processors);
203 if ((relabundfile == "") && (sharedfile == "")) {
204 //is there are current file available for either of these?
205 //give priority to shared, then relabund
206 sharedfile = m->getSharedFile();
207 if (sharedfile != "") { inputFileName = sharedfile; m->mothurOut("Using " + sharedfile + " as input file for the shared parameter."); m->mothurOutEndLine(); }
209 relabundfile = m->getRelAbundFile();
210 if (relabundfile != "") { inputFileName = relabundfile; m->mothurOut("Using " + relabundfile + " as input file for the relabund parameter."); m->mothurOutEndLine(); }
212 m->mothurOut("No valid current files. You must provide a shared or relabund."); m->mothurOutEndLine();
219 if ((designfile == "") && (treefile == "")) {
220 treefile = m->getTreeFile();
221 if (treefile != "") { m->mothurOut("Using " + treefile + " as input file for the tree parameter."); m->mothurOutEndLine(); }
223 designfile = m->getDesignFile();
224 if (designfile != "") { m->mothurOut("Using " + designfile + " as input file for the design parameter."); m->mothurOutEndLine(); }
226 m->mothurOut("[ERROR]: You must provide either a tree or design file."); m->mothurOutEndLine(); abort = true;
231 if ((relabundfile != "") && (sharedfile != "")) { m->mothurOut("[ERROR]: You may not use both a shared and relabund file."); m->mothurOutEndLine(); abort = true; }
236 catch(exception& e) {
237 m->errorOut(e, "IndicatorCommand", "IndicatorCommand");
241 //**********************************************************************************************************************
243 int IndicatorCommand::execute(){
246 if (abort == true) { if (calledHelp) { return 0; } return 2; }
248 cout.setf(ios::fixed, ios::floatfield); cout.setf(ios::showpoint);
250 int start = time(NULL);
252 //read designfile if given and set up groups for read of sharedfiles
253 if (designfile != "") {
254 designMap = new GroupMap(designfile);
255 designMap->readDesignMap();
257 //fill Groups - checks for "all" and for any typo groups
259 vector<string> nameGroups = designMap->getNamesOfGroups();
260 util.setGroups(Groups, nameGroups);
261 designMap->setNamesOfGroups(nameGroups);
263 vector<string> namesSeqs = designMap->getNamesSeqs(Groups);
264 m->setGroups(namesSeqs);
267 /***************************************************/
268 // use smart distancing to get right sharedRabund //
269 /***************************************************/
270 if (sharedfile != "") {
272 if (m->control_pressed) { if (designfile != "") { delete designMap; } for (int i = 0; i < lookup.size(); i++) { delete lookup[i]; } return 0; }
273 if (lookup[0] == NULL) { m->mothurOut("[ERROR] reading shared file."); m->mothurOutEndLine(); return 0; }
276 if (m->control_pressed) { if (designfile != "") { delete designMap; } for (int i = 0; i < lookupFloat.size(); i++) { delete lookupFloat[i]; } return 0; }
277 if (lookupFloat[0] == NULL) { m->mothurOut("[ERROR] reading relabund file."); m->mothurOutEndLine(); return 0; }
280 //reset groups if needed
281 if (designfile != "") { m->setGroups(Groups); }
283 /***************************************************/
284 // reading tree info //
285 /***************************************************/
286 if (treefile != "") {
287 string groupfile = "";
288 m->setTreeFile(treefile);
289 Tree* tree = new Tree(treefile); delete tree; //extracts names from tree to make faked out groupmap
290 treeMap = new TreeMap();
291 bool mismatch = false;
293 for (int i = 0; i < m->Treenames.size(); i++) {
294 //sanity check - is this a group that is not in the sharedfile?
295 if (designfile == "") {
296 if (!(m->inUsersGroups(m->Treenames[i], m->getAllGroups()))) {
297 m->mothurOut("[ERROR]: " + m->Treenames[i] + " is not a group in your shared or relabund file."); m->mothurOutEndLine();
300 treeMap->addSeq(m->Treenames[i], "Group1");
302 vector<string> myGroups; myGroups.push_back(m->Treenames[i]);
303 vector<string> myNames = designMap->getNamesSeqs(myGroups);
305 for(int k = 0; k < myNames.size(); k++) {
306 if (!(m->inUsersGroups(myNames[k], m->getAllGroups()))) {
307 m->mothurOut("[ERROR]: " + myNames[k] + " is not a group in your shared or relabund file."); m->mothurOutEndLine();
311 treeMap->addSeq(m->Treenames[i], "Group1");
315 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; }
317 if (mismatch) { //cleanup and exit
318 if (designfile != "") { delete designMap; }
319 if (sharedfile != "") { for (int i = 0; i < lookup.size(); i++) { delete lookup[i]; } }
320 else { for (int i = 0; i < lookupFloat.size(); i++) { delete lookupFloat[i]; } }
325 read = new ReadNewickTree(treefile);
326 int readOk = read->read(treeMap);
328 if (readOk != 0) { m->mothurOut("Read Terminated."); m->mothurOutEndLine(); delete treeMap; delete read; return 0; }
330 vector<Tree*> T = read->getTrees();
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 for (int i = 0; i < T.size(); i++) { delete T[i]; } delete treeMap; return 0;
341 map<string, string> nameMap;
342 T[0]->assembleTree(nameMap);
344 /***************************************************/
345 // create ouptut tree - respecting pickedGroups //
346 /***************************************************/
347 Tree* outputTree = new Tree(m->getNumGroups(), treeMap);
349 outputTree->getSubTree(T[0], m->getGroups());
350 outputTree->assembleTree(nameMap);
352 //no longer need original tree, we have output tree to use and label
353 for (int i = 0; i < T.size(); i++) { delete T[i]; }
355 if (m->control_pressed) {
356 if (designfile != "") { delete designMap; }
357 if (sharedfile != "") { for (int i = 0; i < lookup.size(); i++) { delete lookup[i]; } }
358 else { for (int i = 0; i < lookupFloat.size(); i++) { delete lookupFloat[i]; } }
359 delete outputTree; delete treeMap; return 0;
362 /***************************************************/
363 // get indicator species values //
364 /***************************************************/
365 GetIndicatorSpecies(outputTree);
366 delete outputTree; delete treeMap;
368 }else { //run with design file only
369 //get indicator species
370 GetIndicatorSpecies();
373 if (designfile != "") { delete designMap; }
374 if (sharedfile != "") { for (int i = 0; i < lookup.size(); i++) { delete lookup[i]; } }
375 else { for (int i = 0; i < lookupFloat.size(); i++) { delete lookupFloat[i]; } }
377 if (m->control_pressed) { for (int i = 0; i < outputNames.size(); i++) { m->mothurRemove(outputNames[i]); } return 0; }
379 //set tree file as new current treefile
380 if (treefile != "") {
382 itTypes = outputTypes.find("tree");
383 if (itTypes != outputTypes.end()) {
384 if ((itTypes->second).size() != 0) { current = (itTypes->second)[0]; m->setTreeFile(current); }
388 m->mothurOutEndLine(); m->mothurOutEndLine(); m->mothurOut("It took " + toString(time(NULL) - start) + " secs to find the indicator species."); m->mothurOutEndLine();
389 m->mothurOutEndLine();
390 m->mothurOut("Output File Names: "); m->mothurOutEndLine();
391 for (int i = 0; i < outputNames.size(); i++) { m->mothurOut(outputNames[i]); m->mothurOutEndLine(); }
392 m->mothurOutEndLine();
396 catch(exception& e) {
397 m->errorOut(e, "IndicatorCommand", "execute");
401 //**********************************************************************************************************************
402 //divide shared or relabund file by groupings in the design file
403 //report all otu values to file
404 int IndicatorCommand::GetIndicatorSpecies(){
406 string thisOutputDir = outputDir;
407 if (outputDir == "") { thisOutputDir += m->hasPath(inputFileName); }
408 string outputFileName = thisOutputDir + m->getRootName(m->getSimpleName(inputFileName)) + getOutputFileNameTag("summary");
409 outputNames.push_back(outputFileName); outputTypes["summary"].push_back(outputFileName);
412 m->openOutputFile(outputFileName, out);
413 out.setf(ios::fixed, ios::floatfield); out.setf(ios::showpoint);
414 m->mothurOutEndLine(); m->mothurOut("Species\tIndicatorValue\tpValue\n");
417 if (sharedfile != "") { numBins = lookup[0]->getNumBins(); }
418 else { numBins = lookupFloat[0]->getNumBins(); }
420 if (m->control_pressed) { out.close(); return 0; }
422 /*****************************************************/
423 //create vectors containing rabund info //
424 /*****************************************************/
426 vector<float> indicatorValues; //size of numBins
427 vector<float> pValues;
428 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.
430 if (sharedfile != "") {
431 vector< vector<SharedRAbundVector*> > groupings;
432 set<string> groupsAlreadyAdded;
433 vector<SharedRAbundVector*> subset;
436 for (int i = 0; i < (designMap->getNamesOfGroups()).size(); i++) {
438 for (int k = 0; k < lookup.size(); k++) {
439 //are you from this grouping?
440 if (designMap->getGroup(lookup[k]->getGroup()) == (designMap->getNamesOfGroups())[i]) {
441 subset.push_back(lookup[k]);
442 groupsAlreadyAdded.insert(lookup[k]->getGroup());
445 if (subset.size() != 0) { groupings.push_back(subset); }
449 if (groupsAlreadyAdded.size() != lookup.size()) { m->mothurOut("[ERROR]: could not make proper groupings."); m->mothurOutEndLine(); }
451 indicatorValues = getValues(groupings, randomGroupingsMap);
453 pValues = getPValues(groupings, randomGroupingsMap, lookup.size(), indicatorValues);
455 vector< vector<SharedRAbundFloatVector*> > groupings;
456 set<string> groupsAlreadyAdded;
457 vector<SharedRAbundFloatVector*> subset;
460 for (int i = 0; i < (designMap->getNamesOfGroups()).size(); i++) {
461 for (int k = 0; k < lookupFloat.size(); k++) {
462 //are you from this grouping?
463 if (designMap->getGroup(lookupFloat[k]->getGroup()) == (designMap->getNamesOfGroups())[i]) {
464 subset.push_back(lookupFloat[k]);
465 groupsAlreadyAdded.insert(lookupFloat[k]->getGroup());
468 if (subset.size() != 0) { groupings.push_back(subset); }
472 if (groupsAlreadyAdded.size() != lookupFloat.size()) { m->mothurOut("[ERROR]: could not make proper groupings."); m->mothurOutEndLine(); }
474 indicatorValues = getValues(groupings, randomGroupingsMap);
476 pValues = getPValues(groupings, randomGroupingsMap, lookupFloat.size(), indicatorValues);
479 if (m->control_pressed) { out.close(); return 0; }
482 /******************************************************/
483 //output indicator values to table form //
484 /*****************************************************/
485 out << "OTU\tIndicator_Value\tpValue" << endl;
486 for (int j = 0; j < indicatorValues.size(); j++) {
488 if (m->control_pressed) { out.close(); return 0; }
490 out << m->currentBinLabels[j] << '\t' << indicatorValues[j] << '\t';
492 if (pValues[j] > (1/(float)iters)) { out << pValues[j] << endl; }
493 else { out << "<" << (1/(float)iters) << endl; }
495 if (pValues[j] <= 0.05) {
496 cout << m->currentBinLabels[j] << '\t' << indicatorValues[j] << '\t';
497 string pValueString = "<" + toString((1/(float)iters));
498 if (pValues[j] > (1/(float)iters)) { pValueString = toString(pValues[j]); cout << pValues[j];}
499 else { cout << "<" << (1/(float)iters); }
500 m->mothurOutJustToLog(m->currentBinLabels[j] + "\t" + toString(indicatorValues[j]) + "\t" + pValueString);
501 m->mothurOutEndLine();
509 catch(exception& e) {
510 m->errorOut(e, "IndicatorCommand", "GetIndicatorSpecies");
514 //**********************************************************************************************************************
515 //traverse tree finding indicator species values for each otu at each node
516 //label node with otu number that has highest indicator value
517 //report all otu values to file
518 int IndicatorCommand::GetIndicatorSpecies(Tree*& T){
521 string thisOutputDir = outputDir;
522 if (outputDir == "") { thisOutputDir += m->hasPath(inputFileName); }
523 string outputFileName = thisOutputDir + m->getRootName(m->getSimpleName(inputFileName)) + getOutputFileNameTag("summary");
524 outputNames.push_back(outputFileName); outputTypes["summary"].push_back(outputFileName);
527 m->openOutputFile(outputFileName, out);
528 out.setf(ios::fixed, ios::floatfield); out.setf(ios::showpoint);
531 if (sharedfile != "") { numBins = lookup[0]->getNumBins(); }
532 else { numBins = lookupFloat[0]->getNumBins(); }
536 for (int i = 0; i < numBins; i++) { out << m->currentBinLabels[i] << "_IndValue" << '\t' << "pValue" << '\t'; }
539 m->mothurOutEndLine(); m->mothurOut("Node\tSpecies\tIndicatorValue\tpValue\n");
541 string treeOutputDir = outputDir;
542 if (outputDir == "") { treeOutputDir += m->hasPath(treefile); }
543 string outputTreeFileName = treeOutputDir + m->getRootName(m->getSimpleName(treefile)) + getOutputFileNameTag("tree");
546 //create a map from tree node index to names of descendants, save time later to know which sharedRabund you need
547 map<int, set<string> > nodeToDescendants;
548 map<int, set<int> > descendantNodes;
549 for (int i = 0; i < T->getNumNodes(); i++) {
550 if (m->control_pressed) { return 0; }
552 nodeToDescendants[i] = getDescendantList(T, i, nodeToDescendants, descendantNodes);
555 //you need the distances to leaf to decide grouping below
556 //this will also set branch lengths if the tree does not include them
557 map<int, float> distToRoot = getDistToRoot(T);
560 for (int i = T->getNumLeaves(); i < T->getNumNodes(); i++) {
561 //cout << endl << i+1 << endl;
562 if (m->control_pressed) { out.close(); return 0; }
564 /*****************************************************/
565 //create vectors containing rabund info //
566 /*****************************************************/
568 vector<float> indicatorValues; //size of numBins
569 vector<float> pValues;
570 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.
572 if (sharedfile != "") {
573 vector< vector<SharedRAbundVector*> > groupings;
575 //get nodes that will be a valid grouping
576 //you are valid if you are not one of my descendants
577 //AND your distToRoot is >= mine
578 //AND you were not added as part of a larger grouping. Largest nodes are added first.
580 set<string> groupsAlreadyAdded;
581 //create a grouping with my grouping
582 vector<SharedRAbundVector*> subset;
584 int doneCount = nodeToDescendants[i].size();
585 for (int k = 0; k < lookup.size(); k++) {
586 //is this descendant of i
587 if ((nodeToDescendants[i].count(lookup[k]->getGroup()) != 0)) {
588 subset.push_back(lookup[k]);
589 groupsAlreadyAdded.insert(lookup[k]->getGroup());
592 if (count == doneCount) { break; } //quit once you get the rabunds for this grouping
594 if (subset.size() != 0) { groupings.push_back(subset); }
597 for (int j = (T->getNumNodes()-1); j >= 0; j--) {
600 if ((descendantNodes[i].count(j) == 0) && (distToRoot[j] >= distToRoot[i])) {
601 vector<SharedRAbundVector*> subset;
603 int doneCount = nodeToDescendants[j].size();
604 for (int k = 0; k < lookup.size(); k++) {
605 //is this descendant of j, and we didn't already add this as part of a larger grouping
606 if ((nodeToDescendants[j].count(lookup[k]->getGroup()) != 0) && (groupsAlreadyAdded.count(lookup[k]->getGroup()) == 0)) {
607 subset.push_back(lookup[k]);
608 groupsAlreadyAdded.insert(lookup[k]->getGroup());
611 if (count == doneCount) { break; } //quit once you get the rabunds for this grouping
614 //if subset.size == 0 then the node was added as part of a larger grouping
615 if (subset.size() != 0) { groupings.push_back(subset); }
619 if (groupsAlreadyAdded.size() != lookup.size()) { m->mothurOut("[ERROR]: could not make proper groupings."); m->mothurOutEndLine(); }
621 indicatorValues = getValues(groupings, randomGroupingsMap);
623 pValues = getPValues(groupings, randomGroupingsMap, lookup.size(), indicatorValues);
625 vector< vector<SharedRAbundFloatVector*> > groupings;
627 //get nodes that will be a valid grouping
628 //you are valid if you are not one of my descendants
629 //AND your distToRoot is >= mine
630 //AND you were not added as part of a larger grouping. Largest nodes are added first.
632 set<string> groupsAlreadyAdded;
633 //create a grouping with my grouping
634 vector<SharedRAbundFloatVector*> subset;
636 int doneCount = nodeToDescendants[i].size();
637 for (int k = 0; k < lookupFloat.size(); k++) {
638 //is this descendant of i
639 if ((nodeToDescendants[i].count(lookupFloat[k]->getGroup()) != 0)) {
640 subset.push_back(lookupFloat[k]);
641 groupsAlreadyAdded.insert(lookupFloat[k]->getGroup());
644 if (count == doneCount) { break; } //quit once you get the rabunds for this grouping
646 if (subset.size() != 0) { groupings.push_back(subset); }
648 for (int j = (T->getNumNodes()-1); j >= 0; j--) {
649 if ((descendantNodes[i].count(j) == 0) && (distToRoot[j] >= distToRoot[i])) {
650 vector<SharedRAbundFloatVector*> subset;
652 int doneCount = nodeToDescendants[j].size();
653 for (int k = 0; k < lookupFloat.size(); k++) {
654 //is this descendant of j, and we didn't already add this as part of a larger grouping
655 if ((nodeToDescendants[j].count(lookupFloat[k]->getGroup()) != 0) && (groupsAlreadyAdded.count(lookupFloat[k]->getGroup()) == 0)) {
656 subset.push_back(lookupFloat[k]);
657 groupsAlreadyAdded.insert(lookupFloat[k]->getGroup());
660 if (count == doneCount) { break; } //quit once you get the rabunds for this grouping
663 //if subset.size == 0 then the node was added as part of a larger grouping
664 if (subset.size() != 0) { groupings.push_back(subset); }
668 if (groupsAlreadyAdded.size() != lookupFloat.size()) { m->mothurOut("[ERROR]: could not make proper groupings."); m->mothurOutEndLine(); }
670 indicatorValues = getValues(groupings, randomGroupingsMap);
672 pValues = getPValues(groupings, randomGroupingsMap, lookupFloat.size(), indicatorValues);
675 if (m->control_pressed) { out.close(); return 0; }
678 /******************************************************/
679 //output indicator values to table form + label tree //
680 /*****************************************************/
681 out << (i+1) << '\t';
682 for (int j = 0; j < indicatorValues.size(); j++) {
684 if (m->control_pressed) { out.close(); return 0; }
686 if (pValues[j] < (1/(float)iters)) {
687 out << indicatorValues[j] << '\t' << '<' << (1/(float)iters) << '\t';
689 out << indicatorValues[j] << '\t' << pValues[j] << '\t';
692 if (pValues[j] <= 0.05) {
693 cout << i+1 << '\t' << m->currentBinLabels[j] << '\t' << indicatorValues[j] << '\t';
694 string pValueString = "<" + toString((1/(float)iters));
695 if (pValues[j] > (1/(float)iters)) { pValueString = toString(pValues[j]); cout << pValues[j];}
696 else { cout << "<" << (1/(float)iters); }
697 m->mothurOutJustToLog(toString(i) + "\t" + m->currentBinLabels[j] + "\t" + toString(indicatorValues[j]) + "\t" + pValueString);
698 m->mothurOutEndLine();
703 T->tree[i].setLabel((i+1));
708 m->openOutputFile(outputTreeFileName, outTree);
709 outputNames.push_back(outputTreeFileName); outputTypes["tree"].push_back(outputTreeFileName);
711 T->print(outTree, "both");
716 catch(exception& e) {
717 m->errorOut(e, "IndicatorCommand", "GetIndicatorSpecies");
721 //**********************************************************************************************************************
722 vector<float> IndicatorCommand::getValues(vector< vector<SharedRAbundFloatVector*> >& groupings, map< vector<int>, vector<int> > groupingsMap){
724 vector<float> values;
725 map< vector<int>, vector<int> >::iterator it;
728 for (int i = 0; i < groupings[0][0]->getNumBins(); i++) {
730 if (m->control_pressed) { return values; }
733 float AijDenominator = 0.0;
736 //get overall abundance of each grouping
737 for (int j = 0; j < groupings.size(); j++) {
739 float totalAbund = 0;
741 for (int k = 0; k < groupings[j].size(); k++) {
742 vector<int> temp; temp.push_back(j); temp.push_back(k);
743 it = groupingsMap.find(temp);
745 if (it == groupingsMap.end()) { //this one didnt get moved
746 totalAbund += groupings[j][k]->getAbundance(i);
747 if (groupings[j][k]->getAbundance(i) != 0.0) { numNotZero++; }
749 totalAbund += groupings[(it->second)[0]][(it->second)[1]]->getAbundance(i);
750 if (groupings[(it->second)[0]][(it->second)[1]]->getAbundance(i) != 0.0) { numNotZero++; }
756 float Aij = (totalAbund / (float) groupings[j].size());
757 terms.push_back(Aij);
759 //percentage of sites represented
760 Bij.push_back(numNotZero / (float) groupings[j].size());
762 AijDenominator += Aij;
765 float maxIndVal = 0.0;
766 for (int j = 0; j < terms.size(); j++) {
767 float thisAij = (terms[j] / AijDenominator); //relative abundance
768 float thisValue = thisAij * Bij[j] * 100.0;
771 if (thisValue > maxIndVal) { maxIndVal = thisValue; }
774 values.push_back(maxIndVal);
779 catch(exception& e) {
780 m->errorOut(e, "IndicatorCommand", "getValues");
784 //**********************************************************************************************************************
785 //same as above, just data type difference
786 vector<float> IndicatorCommand::getValues(vector< vector<SharedRAbundVector*> >& groupings, map< vector<int>, vector<int> > groupingsMap){
788 vector<float> values;
790 /*for (int j = 0; j < groupings.size(); j++) {
791 cout << "grouping " << j << endl;
792 for (int k = 0; k < groupings[j].size(); k++) {
793 cout << groupings[j][k]->getGroup() << endl;
796 map< vector<int>, vector<int> >::iterator it;
799 for (int i = 0; i < groupings[0][0]->getNumBins(); i++) {
801 float AijDenominator = 0.0;
803 //get overall abundance of each grouping
804 for (int j = 0; j < groupings.size(); j++) {
806 int totalAbund = 0.0;
808 for (int k = 0; k < groupings[j].size(); k++) {
809 vector<int> temp; temp.push_back(j); temp.push_back(k);
810 it = groupingsMap.find(temp);
812 if (it == groupingsMap.end()) { //this one didnt get moved
813 totalAbund += groupings[j][k]->getAbundance(i);
814 if (groupings[j][k]->getAbundance(i) != 0.0) { numNotZero++; }
816 totalAbund += groupings[(it->second)[0]][(it->second)[1]]->getAbundance(i);
817 if (groupings[(it->second)[0]][(it->second)[1]]->getAbundance(i) != 0.0) { numNotZero++; }
823 float Aij = (totalAbund / (float) groupings[j].size());
824 terms.push_back(Aij);
826 //percentage of sites represented
827 Bij.push_back(numNotZero / (float) groupings[j].size());
829 AijDenominator += Aij;
832 float maxIndVal = 0.0;
833 for (int j = 0; j < terms.size(); j++) {
834 float thisAij = (terms[j] / AijDenominator); //relative abundance
835 float thisValue = thisAij * Bij[j] * 100.0;
838 if (thisValue > maxIndVal) { maxIndVal = thisValue; }
841 values.push_back(maxIndVal);
846 catch(exception& e) {
847 m->errorOut(e, "IndicatorCommand", "getValues");
851 //**********************************************************************************************************************
852 //you need the distances to root to decide groupings
853 //this will also set branch lengths if the tree does not include them
854 map<int, float> IndicatorCommand::getDistToRoot(Tree*& T){
856 map<int, float> dists;
858 bool hasBranchLengths = false;
859 for (int i = 0; i < T->getNumNodes(); i++) {
860 if (T->tree[i].getBranchLength() > 0.0) { hasBranchLengths = true; break; }
863 //set branchlengths if needed
864 if (!hasBranchLengths) {
865 for (int i = 0; i < T->getNumNodes(); i++) {
867 int lc = T->tree[i].getLChild();
868 int rc = T->tree[i].getRChild();
870 if (lc == -1) { // you are a leaf
871 //if you are a leaf set you priliminary length to 1.0, this may adjust later
872 T->tree[i].setBranchLength(1.0);
874 }else{ // you are an internal node
875 //look at your children's length to leaf
876 float ldist = dists[lc];
877 float rdist = dists[rc];
879 float greater = ldist;
880 if (rdist > greater) { greater = rdist; dists[i] = ldist + 1.0;}
881 else { dists[i] = rdist + 1.0; }
884 //branch length = difference + 1
885 T->tree[lc].setBranchLength((abs(ldist-greater) + 1.0));
886 T->tree[rc].setBranchLength((abs(rdist-greater) + 1.0));
893 for (int i = 0; i < T->getNumNodes(); i++) {
898 while(T->tree[index].getParent() != -1){
899 if (T->tree[index].getBranchLength() != -1) {
900 sum += abs(T->tree[index].getBranchLength());
902 index = T->tree[index].getParent();
910 catch(exception& e) {
911 m->errorOut(e, "IndicatorCommand", "getLengthToLeaf");
915 //**********************************************************************************************************************
916 set<string> IndicatorCommand::getDescendantList(Tree*& T, int i, map<int, set<string> > descendants, map<int, set<int> >& nodes){
920 set<string>::iterator it;
922 int lc = T->tree[i].getLChild();
923 int rc = T->tree[i].getRChild();
925 if (lc == -1) { //you are a leaf your only descendant is yourself
926 set<int> temp; temp.insert(i);
929 if (designfile == "") {
930 names.insert(T->tree[i].getName());
932 vector<string> myGroup; myGroup.push_back(T->tree[i].getName());
933 vector<string> myReps = designMap->getNamesSeqs(myGroup);
934 for (int k = 0; k < myReps.size(); k++) {
935 names.insert(myReps[k]);
939 }else{ //your descedants are the combination of your childrens descendants
940 names = descendants[lc];
941 nodes[i] = nodes[lc];
942 for (it = descendants[rc].begin(); it != descendants[rc].end(); it++) {
945 for (set<int>::iterator itNum = nodes[rc].begin(); itNum != nodes[rc].end(); itNum++) {
946 nodes[i].insert(*itNum);
948 //you are your own descendant
954 catch(exception& e) {
955 m->errorOut(e, "IndicatorCommand", "getDescendantList");
959 //**********************************************************************************************************************
960 int IndicatorCommand::getShared(){
962 InputData* input = new InputData(sharedfile, "sharedfile");
963 lookup = input->getSharedRAbundVectors();
964 string lastLabel = lookup[0]->getLabel();
966 if (label == "") { label = lastLabel; delete input; return 0; }
968 //if the users enters label "0.06" and there is no "0.06" in their file use the next lowest label.
969 set<string> labels; labels.insert(label);
970 set<string> processedLabels;
971 set<string> userLabels = labels;
973 //as long as you are not at the end of the file or done wih the lines you want
974 while((lookup[0] != NULL) && (userLabels.size() != 0)) {
975 if (m->control_pressed) { delete input; return 0; }
977 if(labels.count(lookup[0]->getLabel()) == 1){
978 processedLabels.insert(lookup[0]->getLabel());
979 userLabels.erase(lookup[0]->getLabel());
983 if ((m->anyLabelsToProcess(lookup[0]->getLabel(), userLabels, "") == true) && (processedLabels.count(lastLabel) != 1)) {
984 string saveLabel = lookup[0]->getLabel();
986 for (int i = 0; i < lookup.size(); i++) { delete lookup[i]; }
987 lookup = input->getSharedRAbundVectors(lastLabel);
989 processedLabels.insert(lookup[0]->getLabel());
990 userLabels.erase(lookup[0]->getLabel());
992 //restore real lastlabel to save below
993 lookup[0]->setLabel(saveLabel);
997 lastLabel = lookup[0]->getLabel();
999 //get next line to process
1000 //prevent memory leak
1001 for (int i = 0; i < lookup.size(); i++) { delete lookup[i]; }
1002 lookup = input->getSharedRAbundVectors();
1006 if (m->control_pressed) { delete input; return 0; }
1008 //output error messages about any remaining user labels
1009 set<string>::iterator it;
1010 bool needToRun = false;
1011 for (it = userLabels.begin(); it != userLabels.end(); it++) {
1012 m->mothurOut("Your file does not include the label " + *it);
1013 if (processedLabels.count(lastLabel) != 1) {
1014 m->mothurOut(". I will use " + lastLabel + "."); m->mothurOutEndLine();
1017 m->mothurOut(". Please refer to " + lastLabel + "."); m->mothurOutEndLine();
1021 //run last label if you need to
1022 if (needToRun == true) {
1023 for (int i = 0; i < lookup.size(); i++) { if (lookup[i] != NULL) { delete lookup[i]; } }
1024 lookup = input->getSharedRAbundVectors(lastLabel);
1030 catch(exception& e) {
1031 m->errorOut(e, "IndicatorCommand", "getShared");
1035 //**********************************************************************************************************************
1036 int IndicatorCommand::getSharedFloat(){
1038 InputData* input = new InputData(relabundfile, "relabund");
1039 lookupFloat = input->getSharedRAbundFloatVectors();
1040 string lastLabel = lookupFloat[0]->getLabel();
1042 if (label == "") { label = lastLabel; delete input; return 0; }
1044 //if the users enters label "0.06" and there is no "0.06" in their file use the next lowest label.
1045 set<string> labels; labels.insert(label);
1046 set<string> processedLabels;
1047 set<string> userLabels = labels;
1049 //as long as you are not at the end of the file or done wih the lines you want
1050 while((lookupFloat[0] != NULL) && (userLabels.size() != 0)) {
1052 if (m->control_pressed) { delete input; return 0; }
1054 if(labels.count(lookupFloat[0]->getLabel()) == 1){
1055 processedLabels.insert(lookupFloat[0]->getLabel());
1056 userLabels.erase(lookupFloat[0]->getLabel());
1060 if ((m->anyLabelsToProcess(lookupFloat[0]->getLabel(), userLabels, "") == true) && (processedLabels.count(lastLabel) != 1)) {
1061 string saveLabel = lookupFloat[0]->getLabel();
1063 for (int i = 0; i < lookupFloat.size(); i++) { delete lookupFloat[i]; }
1064 lookupFloat = input->getSharedRAbundFloatVectors(lastLabel);
1066 processedLabels.insert(lookupFloat[0]->getLabel());
1067 userLabels.erase(lookupFloat[0]->getLabel());
1069 //restore real lastlabel to save below
1070 lookupFloat[0]->setLabel(saveLabel);
1074 lastLabel = lookupFloat[0]->getLabel();
1076 //get next line to process
1077 //prevent memory leak
1078 for (int i = 0; i < lookupFloat.size(); i++) { delete lookupFloat[i]; }
1079 lookupFloat = input->getSharedRAbundFloatVectors();
1083 if (m->control_pressed) { delete input; return 0; }
1085 //output error messages about any remaining user labels
1086 set<string>::iterator it;
1087 bool needToRun = false;
1088 for (it = userLabels.begin(); it != userLabels.end(); it++) {
1089 m->mothurOut("Your file does not include the label " + *it);
1090 if (processedLabels.count(lastLabel) != 1) {
1091 m->mothurOut(". I will use " + lastLabel + "."); m->mothurOutEndLine();
1094 m->mothurOut(". Please refer to " + lastLabel + "."); m->mothurOutEndLine();
1098 //run last label if you need to
1099 if (needToRun == true) {
1100 for (int i = 0; i < lookupFloat.size(); i++) { if (lookupFloat[i] != NULL) { delete lookupFloat[i]; } }
1101 lookupFloat = input->getSharedRAbundFloatVectors(lastLabel);
1107 catch(exception& e) {
1108 m->errorOut(e, "IndicatorCommand", "getShared");
1112 //**********************************************************************************************************************
1113 vector<float> IndicatorCommand::driver(vector< vector<SharedRAbundFloatVector*> >& groupings, map< vector<int>, vector<int> > groupingsMap, int num, vector<float> indicatorValues, int numIters){
1115 vector<float> pvalues;
1116 pvalues.resize(indicatorValues.size(), 0);
1118 for(int i=0;i<numIters;i++){
1119 if (m->control_pressed) { break; }
1120 groupingsMap = randomizeGroupings(groupings, num);
1121 vector<float> randomIndicatorValues = getValues(groupings, groupingsMap);
1123 for (int j = 0; j < indicatorValues.size(); j++) {
1124 if (randomIndicatorValues[j] >= indicatorValues[j]) { pvalues[j]++; }
1130 }catch(exception& e) {
1131 m->errorOut(e, "IndicatorCommand", "driver");
1135 //**********************************************************************************************************************
1136 vector<float> IndicatorCommand::getPValues(vector< vector<SharedRAbundFloatVector*> >& groupings, map< vector<int>, vector<int> > groupingsMap, int num, vector<float> indicatorValues){
1138 vector<float> pvalues;
1140 #if defined (__APPLE__) || (__MACH__) || (linux) || (__linux) || (__linux__) || (__unix__) || (__unix)
1141 if(processors == 1){
1142 pvalues = driver(groupings, groupingsMap, num, indicatorValues, iters);
1143 for (int i = 0; i < pvalues.size(); i++) { pvalues[i] /= (double)iters; }
1146 //divide iters between processors
1147 int numItersPerProcessor = iters / processors;
1149 vector<int> processIDS;
1153 //loop through and create all the processes you want
1154 while (process != processors) {
1158 processIDS.push_back(pid); //create map from line number to pid so you can append files in correct order later
1160 }else if (pid == 0){
1161 pvalues = driver(groupings, groupingsMap, num, indicatorValues, numItersPerProcessor);
1163 //pass pvalues to parent
1165 string tempFile = toString(getpid()) + ".pvalues.temp";
1166 m->openOutputFile(tempFile, out);
1169 for (int i = 0; i < pvalues.size(); i++) {
1170 out << pvalues[i] << '\t';
1178 m->mothurOut("[ERROR]: unable to spawn the necessary processes."); m->mothurOutEndLine();
1179 for (int i = 0; i < processIDS.size(); i++) { kill (processIDS[i], SIGINT); }
1185 //special case for last processor in case it doesn't divide evenly
1186 numItersPerProcessor = iters - ((processors-1) * numItersPerProcessor);
1187 pvalues = driver(groupings, groupingsMap, num, indicatorValues, numItersPerProcessor);
1189 //force parent to wait until all the processes are done
1190 for (int i=0;i<processIDS.size();i++) {
1191 int temp = processIDS[i];
1196 for (int i = 0; i < processIDS.size(); i++) {
1198 string tempFile = toString(processIDS[i]) + ".pvalues.temp";
1199 m->openInputFile(tempFile, in);
1201 ////// to do ///////////
1202 int numTemp; numTemp = 0;
1203 for (int j = 0; j < pvalues.size(); j++) {
1204 in >> numTemp; m->gobble(in);
1205 pvalues[j] += numTemp;
1207 in.close(); m->mothurRemove(tempFile);
1209 for (int i = 0; i < pvalues.size(); i++) { pvalues[i] /= (double)iters; }
1212 pvalues = driver(groupings, groupingsMap, num, indicatorValues, iters);
1213 for (int i = 0; i < pvalues.size(); i++) { pvalues[i] /= (double)iters; }
1218 catch(exception& e) {
1219 m->errorOut(e, "IndicatorCommand", "getPValues");
1224 //**********************************************************************************************************************
1225 //same as above, just data type difference
1226 vector<float> IndicatorCommand::driver(vector< vector<SharedRAbundVector*> >& groupings, map< vector<int>, vector<int> > groupingsMap, int num, vector<float> indicatorValues, int numIters){
1228 vector<float> pvalues;
1229 pvalues.resize(indicatorValues.size(), 0);
1231 for(int i=0;i<numIters;i++){
1232 if (m->control_pressed) { break; }
1233 groupingsMap = randomizeGroupings(groupings, num);
1234 vector<float> randomIndicatorValues = getValues(groupings, groupingsMap);
1236 for (int j = 0; j < indicatorValues.size(); j++) {
1237 if (randomIndicatorValues[j] >= indicatorValues[j]) { pvalues[j]++; }
1243 }catch(exception& e) {
1244 m->errorOut(e, "IndicatorCommand", "driver");
1248 //**********************************************************************************************************************
1249 //same as above, just data type difference
1250 vector<float> IndicatorCommand::getPValues(vector< vector<SharedRAbundVector*> >& groupings, map< vector<int>, vector<int> > groupingsMap, int num, vector<float> indicatorValues){
1252 vector<float> pvalues;
1254 #if defined (__APPLE__) || (__MACH__) || (linux) || (__linux) || (__linux__) || (__unix__) || (__unix)
1255 if(processors == 1){
1256 pvalues = driver(groupings, groupingsMap, num, indicatorValues, iters);
1257 for (int i = 0; i < pvalues.size(); i++) { pvalues[i] /= (double)iters; }
1260 //divide iters between processors
1261 int numItersPerProcessor = iters / processors;
1263 vector<int> processIDS;
1266 //loop through and create all the processes you want
1267 while (process != processors) {
1271 processIDS.push_back(pid); //create map from line number to pid so you can append files in correct order later
1273 }else if (pid == 0){
1274 pvalues = driver(groupings, groupingsMap, num, indicatorValues, numItersPerProcessor);
1276 //pass pvalues to parent
1278 string tempFile = toString(getpid()) + ".pvalues.temp";
1279 m->openOutputFile(tempFile, out);
1282 for (int i = 0; i < pvalues.size(); i++) {
1283 out << pvalues[i] << '\t';
1291 m->mothurOut("[ERROR]: unable to spawn the necessary processes."); m->mothurOutEndLine();
1292 for (int i = 0; i < processIDS.size(); i++) { kill (processIDS[i], SIGINT); }
1298 //special case for last processor in case it doesn't divide evenly
1299 numItersPerProcessor = iters - ((processors-1) * numItersPerProcessor);
1300 pvalues = driver(groupings, groupingsMap, num, indicatorValues, numItersPerProcessor);
1302 //force parent to wait until all the processes are done
1303 for (int i=0;i<processIDS.size();i++) {
1304 int temp = processIDS[i];
1309 for (int i = 0; i < processIDS.size(); i++) {
1311 string tempFile = toString(processIDS[i]) + ".pvalues.temp";
1312 m->openInputFile(tempFile, in);
1314 ////// to do ///////////
1315 int numTemp; numTemp = 0;
1316 for (int j = 0; j < pvalues.size(); j++) {
1317 in >> numTemp; m->gobble(in);
1318 pvalues[j] += numTemp;
1320 in.close(); m->mothurRemove(tempFile);
1322 for (int i = 0; i < pvalues.size(); i++) { pvalues[i] /= (double)iters; }
1325 pvalues = driver(groupings, groupingsMap, num, indicatorValues, iters);
1326 for (int i = 0; i < pvalues.size(); i++) { pvalues[i] /= (double)iters; }
1331 catch(exception& e) {
1332 m->errorOut(e, "IndicatorCommand", "getPValues");
1336 //**********************************************************************************************************************
1337 //swap groups between groupings, in essence randomizing the second column of the design file
1338 map< vector<int>, vector<int> > IndicatorCommand::randomizeGroupings(vector< vector<SharedRAbundVector*> >& groupings, int numLookupGroups){
1341 map< vector<int>, vector<int> > randomGroupings;
1343 for (int i = 0; i < numLookupGroups; i++) {
1344 if (m->control_pressed) {break;}
1346 //get random groups to swap to switch with
1347 //generate random int between 0 and groupings.size()-1
1348 int z = m->getRandomIndex(groupings.size()-1);
1349 int x = m->getRandomIndex(groupings.size()-1);
1350 int a = m->getRandomIndex(groupings[z].size()-1);
1351 int b = m->getRandomIndex(groupings[x].size()-1);
1352 //cout << i << '\t' << z << '\t' << x << '\t' << a << '\t' << b << endl;
1353 //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; }
1358 from.push_back(z); from.push_back(a);
1359 to.push_back(x); to.push_back(b);
1361 randomGroupings[from] = to;
1363 //cout << "done" << endl;
1364 return randomGroupings;
1366 catch(exception& e) {
1367 m->errorOut(e, "IndicatorCommand", "randomizeGroupings");
1371 //**********************************************************************************************************************
1372 //swap groups between groupings, in essence randomizing the second column of the design file
1373 map< vector<int>, vector<int> > IndicatorCommand::randomizeGroupings(vector< vector<SharedRAbundFloatVector*> >& groupings, int numLookupGroups){
1376 map< vector<int>, vector<int> > randomGroupings;
1378 for (int i = 0; i < numLookupGroups; i++) {
1380 //get random groups to swap to switch with
1381 //generate random int between 0 and groupings.size()-1
1382 int z = m->getRandomIndex(groupings.size()-1);
1383 int x = m->getRandomIndex(groupings.size()-1);
1384 int a = m->getRandomIndex(groupings[z].size()-1);
1385 int b = m->getRandomIndex(groupings[x].size()-1);
1386 //cout << i << '\t' << z << '\t' << x << '\t' << a << '\t' << b << endl;
1391 from.push_back(z); from.push_back(a);
1392 to.push_back(x); to.push_back(b);
1394 randomGroupings[from] = to;
1397 return randomGroupings;
1399 catch(exception& e) {
1400 m->errorOut(e, "IndicatorCommand", "randomizeGroupings");
1405 /*****************************************************************/