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 ct = new CountTable();
291 bool mismatch = false;
294 map<string, string> groupMap;
296 for (int i = 0; i < m->Treenames.size(); i++) {
297 nameMap.insert(m->Treenames[i]);
298 //sanity check - is this a group that is not in the sharedfile?
299 if (designfile == "") {
300 if (i == 0) { gps.insert("Group1"); }
301 if (!(m->inUsersGroups(m->Treenames[i], m->getAllGroups()))) {
302 m->mothurOut("[ERROR]: " + m->Treenames[i] + " is not a group in your shared or relabund file."); m->mothurOutEndLine();
305 groupMap[m->Treenames[i]] = "Group1";
307 vector<string> myGroups; myGroups.push_back(m->Treenames[i]);
308 vector<string> myNames = designMap->getNamesSeqs(myGroups);
310 for(int k = 0; k < myNames.size(); k++) {
311 if (!(m->inUsersGroups(myNames[k], m->getAllGroups()))) {
312 m->mothurOut("[ERROR]: " + myNames[k] + " is not a group in your shared or relabund file."); m->mothurOutEndLine();
316 groupMap[m->Treenames[i]] = "Group1";
319 ct->createTable(nameMap, groupMap, gps);
321 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; }
323 if (mismatch) { //cleanup and exit
324 if (designfile != "") { delete designMap; }
325 if (sharedfile != "") { for (int i = 0; i < lookup.size(); i++) { delete lookup[i]; } }
326 else { for (int i = 0; i < lookupFloat.size(); i++) { delete lookupFloat[i]; } }
331 read = new ReadNewickTree(treefile);
332 int readOk = read->read(ct);
334 if (readOk != 0) { m->mothurOut("Read Terminated."); m->mothurOutEndLine(); delete ct; delete read; return 0; }
336 vector<Tree*> T = read->getTrees();
340 if (m->control_pressed) {
341 if (designfile != "") { delete designMap; }
342 if (sharedfile != "") { for (int i = 0; i < lookup.size(); i++) { delete lookup[i]; } }
343 else { for (int i = 0; i < lookupFloat.size(); i++) { delete lookupFloat[i]; } }
344 for (int i = 0; i < T.size(); i++) { delete T[i]; } delete ct; return 0;
347 T[0]->assembleTree();
349 /***************************************************/
350 // create ouptut tree - respecting pickedGroups //
351 /***************************************************/
352 Tree* outputTree = new Tree(m->getNumGroups(), ct);
354 outputTree->getSubTree(T[0], m->getGroups());
355 outputTree->assembleTree();
357 //no longer need original tree, we have output tree to use and label
358 for (int i = 0; i < T.size(); i++) { delete T[i]; }
360 if (m->control_pressed) {
361 if (designfile != "") { delete designMap; }
362 if (sharedfile != "") { for (int i = 0; i < lookup.size(); i++) { delete lookup[i]; } }
363 else { for (int i = 0; i < lookupFloat.size(); i++) { delete lookupFloat[i]; } }
364 delete outputTree; delete ct; return 0;
367 /***************************************************/
368 // get indicator species values //
369 /***************************************************/
370 GetIndicatorSpecies(outputTree);
371 delete outputTree; delete ct;
373 }else { //run with design file only
374 //get indicator species
375 GetIndicatorSpecies();
378 if (designfile != "") { delete designMap; }
379 if (sharedfile != "") { for (int i = 0; i < lookup.size(); i++) { delete lookup[i]; } }
380 else { for (int i = 0; i < lookupFloat.size(); i++) { delete lookupFloat[i]; } }
382 if (m->control_pressed) { for (int i = 0; i < outputNames.size(); i++) { m->mothurRemove(outputNames[i]); } return 0; }
384 //set tree file as new current treefile
385 if (treefile != "") {
387 itTypes = outputTypes.find("tree");
388 if (itTypes != outputTypes.end()) {
389 if ((itTypes->second).size() != 0) { current = (itTypes->second)[0]; m->setTreeFile(current); }
393 m->mothurOutEndLine(); m->mothurOutEndLine(); m->mothurOut("It took " + toString(time(NULL) - start) + " secs to find the indicator species."); m->mothurOutEndLine();
394 m->mothurOutEndLine();
395 m->mothurOut("Output File Names: "); m->mothurOutEndLine();
396 for (int i = 0; i < outputNames.size(); i++) { m->mothurOut(outputNames[i]); m->mothurOutEndLine(); }
397 m->mothurOutEndLine();
401 catch(exception& e) {
402 m->errorOut(e, "IndicatorCommand", "execute");
406 //**********************************************************************************************************************
407 //divide shared or relabund file by groupings in the design file
408 //report all otu values to file
409 int IndicatorCommand::GetIndicatorSpecies(){
411 string thisOutputDir = outputDir;
412 if (outputDir == "") { thisOutputDir += m->hasPath(inputFileName); }
413 string outputFileName = thisOutputDir + m->getRootName(m->getSimpleName(inputFileName)) + getOutputFileNameTag("summary");
414 outputNames.push_back(outputFileName); outputTypes["summary"].push_back(outputFileName);
417 m->openOutputFile(outputFileName, out);
418 out.setf(ios::fixed, ios::floatfield); out.setf(ios::showpoint);
419 m->mothurOutEndLine(); m->mothurOut("Species\tIndicatorValue\tpValue\n");
422 if (sharedfile != "") { numBins = lookup[0]->getNumBins(); }
423 else { numBins = lookupFloat[0]->getNumBins(); }
425 if (m->control_pressed) { out.close(); return 0; }
427 /*****************************************************/
428 //create vectors containing rabund info //
429 /*****************************************************/
431 vector<float> indicatorValues; //size of numBins
432 vector<float> pValues;
433 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.
435 if (sharedfile != "") {
436 vector< vector<SharedRAbundVector*> > groupings;
437 set<string> groupsAlreadyAdded;
438 vector<SharedRAbundVector*> subset;
441 for (int i = 0; i < (designMap->getNamesOfGroups()).size(); i++) {
443 for (int k = 0; k < lookup.size(); k++) {
444 //are you from this grouping?
445 if (designMap->getGroup(lookup[k]->getGroup()) == (designMap->getNamesOfGroups())[i]) {
446 subset.push_back(lookup[k]);
447 groupsAlreadyAdded.insert(lookup[k]->getGroup());
450 if (subset.size() != 0) { groupings.push_back(subset); }
454 if (groupsAlreadyAdded.size() != lookup.size()) { m->mothurOut("[ERROR]: could not make proper groupings."); m->mothurOutEndLine(); }
456 indicatorValues = getValues(groupings, randomGroupingsMap);
458 pValues = getPValues(groupings, randomGroupingsMap, lookup.size(), indicatorValues);
460 vector< vector<SharedRAbundFloatVector*> > groupings;
461 set<string> groupsAlreadyAdded;
462 vector<SharedRAbundFloatVector*> subset;
465 for (int i = 0; i < (designMap->getNamesOfGroups()).size(); i++) {
466 for (int k = 0; k < lookupFloat.size(); k++) {
467 //are you from this grouping?
468 if (designMap->getGroup(lookupFloat[k]->getGroup()) == (designMap->getNamesOfGroups())[i]) {
469 subset.push_back(lookupFloat[k]);
470 groupsAlreadyAdded.insert(lookupFloat[k]->getGroup());
473 if (subset.size() != 0) { groupings.push_back(subset); }
477 if (groupsAlreadyAdded.size() != lookupFloat.size()) { m->mothurOut("[ERROR]: could not make proper groupings."); m->mothurOutEndLine(); }
479 indicatorValues = getValues(groupings, randomGroupingsMap);
481 pValues = getPValues(groupings, randomGroupingsMap, lookupFloat.size(), indicatorValues);
484 if (m->control_pressed) { out.close(); return 0; }
487 /******************************************************/
488 //output indicator values to table form //
489 /*****************************************************/
490 out << "OTU\tIndicator_Value\tpValue" << endl;
491 for (int j = 0; j < indicatorValues.size(); j++) {
493 if (m->control_pressed) { out.close(); return 0; }
495 out << m->currentBinLabels[j] << '\t' << indicatorValues[j] << '\t';
497 if (pValues[j] > (1/(float)iters)) { out << pValues[j] << endl; }
498 else { out << "<" << (1/(float)iters) << endl; }
500 if (pValues[j] <= 0.05) {
501 cout << m->currentBinLabels[j] << '\t' << indicatorValues[j] << '\t';
502 string pValueString = "<" + toString((1/(float)iters));
503 if (pValues[j] > (1/(float)iters)) { pValueString = toString(pValues[j]); cout << pValues[j];}
504 else { cout << "<" << (1/(float)iters); }
505 m->mothurOutJustToLog(m->currentBinLabels[j] + "\t" + toString(indicatorValues[j]) + "\t" + pValueString);
506 m->mothurOutEndLine();
514 catch(exception& e) {
515 m->errorOut(e, "IndicatorCommand", "GetIndicatorSpecies");
519 //**********************************************************************************************************************
520 //traverse tree finding indicator species values for each otu at each node
521 //label node with otu number that has highest indicator value
522 //report all otu values to file
523 int IndicatorCommand::GetIndicatorSpecies(Tree*& T){
526 string thisOutputDir = outputDir;
527 if (outputDir == "") { thisOutputDir += m->hasPath(inputFileName); }
528 string outputFileName = thisOutputDir + m->getRootName(m->getSimpleName(inputFileName)) + getOutputFileNameTag("summary");
529 outputNames.push_back(outputFileName); outputTypes["summary"].push_back(outputFileName);
532 m->openOutputFile(outputFileName, out);
533 out.setf(ios::fixed, ios::floatfield); out.setf(ios::showpoint);
536 if (sharedfile != "") { numBins = lookup[0]->getNumBins(); }
537 else { numBins = lookupFloat[0]->getNumBins(); }
541 for (int i = 0; i < numBins; i++) { out << m->currentBinLabels[i] << "_IndValue" << '\t' << "pValue" << '\t'; }
544 m->mothurOutEndLine(); m->mothurOut("Node\tSpecies\tIndicatorValue\tpValue\n");
546 string treeOutputDir = outputDir;
547 if (outputDir == "") { treeOutputDir += m->hasPath(treefile); }
548 string outputTreeFileName = treeOutputDir + m->getRootName(m->getSimpleName(treefile)) + getOutputFileNameTag("tree");
551 //create a map from tree node index to names of descendants, save time later to know which sharedRabund you need
552 map<int, set<string> > nodeToDescendants;
553 map<int, set<int> > descendantNodes;
554 for (int i = 0; i < T->getNumNodes(); i++) {
555 if (m->control_pressed) { return 0; }
557 nodeToDescendants[i] = getDescendantList(T, i, nodeToDescendants, descendantNodes);
560 //you need the distances to leaf to decide grouping below
561 //this will also set branch lengths if the tree does not include them
562 map<int, float> distToRoot = getDistToRoot(T);
565 for (int i = T->getNumLeaves(); i < T->getNumNodes(); i++) {
566 //cout << endl << i+1 << endl;
567 if (m->control_pressed) { out.close(); return 0; }
569 /*****************************************************/
570 //create vectors containing rabund info //
571 /*****************************************************/
573 vector<float> indicatorValues; //size of numBins
574 vector<float> pValues;
575 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.
577 if (sharedfile != "") {
578 vector< vector<SharedRAbundVector*> > groupings;
580 //get nodes that will be a valid grouping
581 //you are valid if you are not one of my descendants
582 //AND your distToRoot is >= mine
583 //AND you were not added as part of a larger grouping. Largest nodes are added first.
585 set<string> groupsAlreadyAdded;
586 //create a grouping with my grouping
587 vector<SharedRAbundVector*> subset;
589 int doneCount = nodeToDescendants[i].size();
590 for (int k = 0; k < lookup.size(); k++) {
591 //is this descendant of i
592 if ((nodeToDescendants[i].count(lookup[k]->getGroup()) != 0)) {
593 subset.push_back(lookup[k]);
594 groupsAlreadyAdded.insert(lookup[k]->getGroup());
597 if (count == doneCount) { break; } //quit once you get the rabunds for this grouping
599 if (subset.size() != 0) { groupings.push_back(subset); }
602 for (int j = (T->getNumNodes()-1); j >= 0; j--) {
605 if ((descendantNodes[i].count(j) == 0) && (distToRoot[j] >= distToRoot[i])) {
606 vector<SharedRAbundVector*> subset;
608 int doneCount = nodeToDescendants[j].size();
609 for (int k = 0; k < lookup.size(); k++) {
610 //is this descendant of j, and we didn't already add this as part of a larger grouping
611 if ((nodeToDescendants[j].count(lookup[k]->getGroup()) != 0) && (groupsAlreadyAdded.count(lookup[k]->getGroup()) == 0)) {
612 subset.push_back(lookup[k]);
613 groupsAlreadyAdded.insert(lookup[k]->getGroup());
616 if (count == doneCount) { break; } //quit once you get the rabunds for this grouping
619 //if subset.size == 0 then the node was added as part of a larger grouping
620 if (subset.size() != 0) { groupings.push_back(subset); }
624 if (groupsAlreadyAdded.size() != lookup.size()) { m->mothurOut("[ERROR]: could not make proper groupings."); m->mothurOutEndLine(); }
626 indicatorValues = getValues(groupings, randomGroupingsMap);
628 pValues = getPValues(groupings, randomGroupingsMap, lookup.size(), indicatorValues);
630 vector< vector<SharedRAbundFloatVector*> > groupings;
632 //get nodes that will be a valid grouping
633 //you are valid if you are not one of my descendants
634 //AND your distToRoot is >= mine
635 //AND you were not added as part of a larger grouping. Largest nodes are added first.
637 set<string> groupsAlreadyAdded;
638 //create a grouping with my grouping
639 vector<SharedRAbundFloatVector*> subset;
641 int doneCount = nodeToDescendants[i].size();
642 for (int k = 0; k < lookupFloat.size(); k++) {
643 //is this descendant of i
644 if ((nodeToDescendants[i].count(lookupFloat[k]->getGroup()) != 0)) {
645 subset.push_back(lookupFloat[k]);
646 groupsAlreadyAdded.insert(lookupFloat[k]->getGroup());
649 if (count == doneCount) { break; } //quit once you get the rabunds for this grouping
651 if (subset.size() != 0) { groupings.push_back(subset); }
653 for (int j = (T->getNumNodes()-1); j >= 0; j--) {
654 if ((descendantNodes[i].count(j) == 0) && (distToRoot[j] >= distToRoot[i])) {
655 vector<SharedRAbundFloatVector*> subset;
657 int doneCount = nodeToDescendants[j].size();
658 for (int k = 0; k < lookupFloat.size(); k++) {
659 //is this descendant of j, and we didn't already add this as part of a larger grouping
660 if ((nodeToDescendants[j].count(lookupFloat[k]->getGroup()) != 0) && (groupsAlreadyAdded.count(lookupFloat[k]->getGroup()) == 0)) {
661 subset.push_back(lookupFloat[k]);
662 groupsAlreadyAdded.insert(lookupFloat[k]->getGroup());
665 if (count == doneCount) { break; } //quit once you get the rabunds for this grouping
668 //if subset.size == 0 then the node was added as part of a larger grouping
669 if (subset.size() != 0) { groupings.push_back(subset); }
673 if (groupsAlreadyAdded.size() != lookupFloat.size()) { m->mothurOut("[ERROR]: could not make proper groupings."); m->mothurOutEndLine(); }
675 indicatorValues = getValues(groupings, randomGroupingsMap);
677 pValues = getPValues(groupings, randomGroupingsMap, lookupFloat.size(), indicatorValues);
680 if (m->control_pressed) { out.close(); return 0; }
683 /******************************************************/
684 //output indicator values to table form + label tree //
685 /*****************************************************/
686 out << (i+1) << '\t';
687 for (int j = 0; j < indicatorValues.size(); j++) {
689 if (m->control_pressed) { out.close(); return 0; }
691 if (pValues[j] < (1/(float)iters)) {
692 out << indicatorValues[j] << '\t' << '<' << (1/(float)iters) << '\t';
694 out << indicatorValues[j] << '\t' << pValues[j] << '\t';
697 if (pValues[j] <= 0.05) {
698 cout << i+1 << '\t' << m->currentBinLabels[j] << '\t' << indicatorValues[j] << '\t';
699 string pValueString = "<" + toString((1/(float)iters));
700 if (pValues[j] > (1/(float)iters)) { pValueString = toString(pValues[j]); cout << pValues[j];}
701 else { cout << "<" << (1/(float)iters); }
702 m->mothurOutJustToLog(toString(i) + "\t" + m->currentBinLabels[j] + "\t" + toString(indicatorValues[j]) + "\t" + pValueString);
703 m->mothurOutEndLine();
708 T->tree[i].setLabel((i+1));
713 m->openOutputFile(outputTreeFileName, outTree);
714 outputNames.push_back(outputTreeFileName); outputTypes["tree"].push_back(outputTreeFileName);
716 T->print(outTree, "both");
721 catch(exception& e) {
722 m->errorOut(e, "IndicatorCommand", "GetIndicatorSpecies");
726 //**********************************************************************************************************************
727 vector<float> IndicatorCommand::getValues(vector< vector<SharedRAbundFloatVector*> >& groupings, map< vector<int>, vector<int> > groupingsMap){
729 vector<float> values;
730 map< vector<int>, vector<int> >::iterator it;
733 for (int i = 0; i < groupings[0][0]->getNumBins(); i++) {
735 if (m->control_pressed) { return values; }
738 float AijDenominator = 0.0;
741 //get overall abundance of each grouping
742 for (int j = 0; j < groupings.size(); j++) {
744 float totalAbund = 0;
746 for (int k = 0; k < groupings[j].size(); k++) {
747 vector<int> temp; temp.push_back(j); temp.push_back(k);
748 it = groupingsMap.find(temp);
750 if (it == groupingsMap.end()) { //this one didnt get moved
751 totalAbund += groupings[j][k]->getAbundance(i);
752 if (groupings[j][k]->getAbundance(i) != 0.0) { numNotZero++; }
754 totalAbund += groupings[(it->second)[0]][(it->second)[1]]->getAbundance(i);
755 if (groupings[(it->second)[0]][(it->second)[1]]->getAbundance(i) != 0.0) { numNotZero++; }
761 float Aij = (totalAbund / (float) groupings[j].size());
762 terms.push_back(Aij);
764 //percentage of sites represented
765 Bij.push_back(numNotZero / (float) groupings[j].size());
767 AijDenominator += Aij;
770 float maxIndVal = 0.0;
771 for (int j = 0; j < terms.size(); j++) {
772 float thisAij = (terms[j] / AijDenominator); //relative abundance
773 float thisValue = thisAij * Bij[j] * 100.0;
776 if (thisValue > maxIndVal) { maxIndVal = thisValue; }
779 values.push_back(maxIndVal);
784 catch(exception& e) {
785 m->errorOut(e, "IndicatorCommand", "getValues");
789 //**********************************************************************************************************************
790 //same as above, just data type difference
791 vector<float> IndicatorCommand::getValues(vector< vector<SharedRAbundVector*> >& groupings, map< vector<int>, vector<int> > groupingsMap){
793 vector<float> values;
795 /*for (int j = 0; j < groupings.size(); j++) {
796 cout << "grouping " << j << endl;
797 for (int k = 0; k < groupings[j].size(); k++) {
798 cout << groupings[j][k]->getGroup() << endl;
801 map< vector<int>, vector<int> >::iterator it;
804 for (int i = 0; i < groupings[0][0]->getNumBins(); i++) {
806 float AijDenominator = 0.0;
808 //get overall abundance of each grouping
809 for (int j = 0; j < groupings.size(); j++) {
811 int totalAbund = 0.0;
813 for (int k = 0; k < groupings[j].size(); k++) {
814 vector<int> temp; temp.push_back(j); temp.push_back(k);
815 it = groupingsMap.find(temp);
817 if (it == groupingsMap.end()) { //this one didnt get moved
818 totalAbund += groupings[j][k]->getAbundance(i);
819 if (groupings[j][k]->getAbundance(i) != 0.0) { numNotZero++; }
821 totalAbund += groupings[(it->second)[0]][(it->second)[1]]->getAbundance(i);
822 if (groupings[(it->second)[0]][(it->second)[1]]->getAbundance(i) != 0.0) { numNotZero++; }
828 float Aij = (totalAbund / (float) groupings[j].size());
829 terms.push_back(Aij);
831 //percentage of sites represented
832 Bij.push_back(numNotZero / (float) groupings[j].size());
834 AijDenominator += Aij;
837 float maxIndVal = 0.0;
838 for (int j = 0; j < terms.size(); j++) {
839 float thisAij = (terms[j] / AijDenominator); //relative abundance
840 float thisValue = thisAij * Bij[j] * 100.0;
843 if (thisValue > maxIndVal) { maxIndVal = thisValue; }
846 values.push_back(maxIndVal);
851 catch(exception& e) {
852 m->errorOut(e, "IndicatorCommand", "getValues");
856 //**********************************************************************************************************************
857 //you need the distances to root to decide groupings
858 //this will also set branch lengths if the tree does not include them
859 map<int, float> IndicatorCommand::getDistToRoot(Tree*& T){
861 map<int, float> dists;
863 bool hasBranchLengths = false;
864 for (int i = 0; i < T->getNumNodes(); i++) {
865 if (T->tree[i].getBranchLength() > 0.0) { hasBranchLengths = true; break; }
868 //set branchlengths if needed
869 if (!hasBranchLengths) {
870 for (int i = 0; i < T->getNumNodes(); i++) {
872 int lc = T->tree[i].getLChild();
873 int rc = T->tree[i].getRChild();
875 if (lc == -1) { // you are a leaf
876 //if you are a leaf set you priliminary length to 1.0, this may adjust later
877 T->tree[i].setBranchLength(1.0);
879 }else{ // you are an internal node
880 //look at your children's length to leaf
881 float ldist = dists[lc];
882 float rdist = dists[rc];
884 float greater = ldist;
885 if (rdist > greater) { greater = rdist; dists[i] = ldist + 1.0;}
886 else { dists[i] = rdist + 1.0; }
889 //branch length = difference + 1
890 T->tree[lc].setBranchLength((abs(ldist-greater) + 1.0));
891 T->tree[rc].setBranchLength((abs(rdist-greater) + 1.0));
898 for (int i = 0; i < T->getNumNodes(); i++) {
903 while(T->tree[index].getParent() != -1){
904 if (T->tree[index].getBranchLength() != -1) {
905 sum += abs(T->tree[index].getBranchLength());
907 index = T->tree[index].getParent();
915 catch(exception& e) {
916 m->errorOut(e, "IndicatorCommand", "getLengthToLeaf");
920 //**********************************************************************************************************************
921 set<string> IndicatorCommand::getDescendantList(Tree*& T, int i, map<int, set<string> > descendants, map<int, set<int> >& nodes){
925 set<string>::iterator it;
927 int lc = T->tree[i].getLChild();
928 int rc = T->tree[i].getRChild();
930 if (lc == -1) { //you are a leaf your only descendant is yourself
931 set<int> temp; temp.insert(i);
934 if (designfile == "") {
935 names.insert(T->tree[i].getName());
937 vector<string> myGroup; myGroup.push_back(T->tree[i].getName());
938 vector<string> myReps = designMap->getNamesSeqs(myGroup);
939 for (int k = 0; k < myReps.size(); k++) {
940 names.insert(myReps[k]);
944 }else{ //your descedants are the combination of your childrens descendants
945 names = descendants[lc];
946 nodes[i] = nodes[lc];
947 for (it = descendants[rc].begin(); it != descendants[rc].end(); it++) {
950 for (set<int>::iterator itNum = nodes[rc].begin(); itNum != nodes[rc].end(); itNum++) {
951 nodes[i].insert(*itNum);
953 //you are your own descendant
959 catch(exception& e) {
960 m->errorOut(e, "IndicatorCommand", "getDescendantList");
964 //**********************************************************************************************************************
965 int IndicatorCommand::getShared(){
967 InputData* input = new InputData(sharedfile, "sharedfile");
968 lookup = input->getSharedRAbundVectors();
969 string lastLabel = lookup[0]->getLabel();
971 if (label == "") { label = lastLabel; delete input; return 0; }
973 //if the users enters label "0.06" and there is no "0.06" in their file use the next lowest label.
974 set<string> labels; labels.insert(label);
975 set<string> processedLabels;
976 set<string> userLabels = labels;
978 //as long as you are not at the end of the file or done wih the lines you want
979 while((lookup[0] != NULL) && (userLabels.size() != 0)) {
980 if (m->control_pressed) { delete input; return 0; }
982 if(labels.count(lookup[0]->getLabel()) == 1){
983 processedLabels.insert(lookup[0]->getLabel());
984 userLabels.erase(lookup[0]->getLabel());
988 if ((m->anyLabelsToProcess(lookup[0]->getLabel(), userLabels, "") == true) && (processedLabels.count(lastLabel) != 1)) {
989 string saveLabel = lookup[0]->getLabel();
991 for (int i = 0; i < lookup.size(); i++) { delete lookup[i]; }
992 lookup = input->getSharedRAbundVectors(lastLabel);
994 processedLabels.insert(lookup[0]->getLabel());
995 userLabels.erase(lookup[0]->getLabel());
997 //restore real lastlabel to save below
998 lookup[0]->setLabel(saveLabel);
1002 lastLabel = lookup[0]->getLabel();
1004 //get next line to process
1005 //prevent memory leak
1006 for (int i = 0; i < lookup.size(); i++) { delete lookup[i]; }
1007 lookup = input->getSharedRAbundVectors();
1011 if (m->control_pressed) { delete input; return 0; }
1013 //output error messages about any remaining user labels
1014 set<string>::iterator it;
1015 bool needToRun = false;
1016 for (it = userLabels.begin(); it != userLabels.end(); it++) {
1017 m->mothurOut("Your file does not include the label " + *it);
1018 if (processedLabels.count(lastLabel) != 1) {
1019 m->mothurOut(". I will use " + lastLabel + "."); m->mothurOutEndLine();
1022 m->mothurOut(". Please refer to " + lastLabel + "."); m->mothurOutEndLine();
1026 //run last label if you need to
1027 if (needToRun == true) {
1028 for (int i = 0; i < lookup.size(); i++) { if (lookup[i] != NULL) { delete lookup[i]; } }
1029 lookup = input->getSharedRAbundVectors(lastLabel);
1035 catch(exception& e) {
1036 m->errorOut(e, "IndicatorCommand", "getShared");
1040 //**********************************************************************************************************************
1041 int IndicatorCommand::getSharedFloat(){
1043 InputData* input = new InputData(relabundfile, "relabund");
1044 lookupFloat = input->getSharedRAbundFloatVectors();
1045 string lastLabel = lookupFloat[0]->getLabel();
1047 if (label == "") { label = lastLabel; delete input; return 0; }
1049 //if the users enters label "0.06" and there is no "0.06" in their file use the next lowest label.
1050 set<string> labels; labels.insert(label);
1051 set<string> processedLabels;
1052 set<string> userLabels = labels;
1054 //as long as you are not at the end of the file or done wih the lines you want
1055 while((lookupFloat[0] != NULL) && (userLabels.size() != 0)) {
1057 if (m->control_pressed) { delete input; return 0; }
1059 if(labels.count(lookupFloat[0]->getLabel()) == 1){
1060 processedLabels.insert(lookupFloat[0]->getLabel());
1061 userLabels.erase(lookupFloat[0]->getLabel());
1065 if ((m->anyLabelsToProcess(lookupFloat[0]->getLabel(), userLabels, "") == true) && (processedLabels.count(lastLabel) != 1)) {
1066 string saveLabel = lookupFloat[0]->getLabel();
1068 for (int i = 0; i < lookupFloat.size(); i++) { delete lookupFloat[i]; }
1069 lookupFloat = input->getSharedRAbundFloatVectors(lastLabel);
1071 processedLabels.insert(lookupFloat[0]->getLabel());
1072 userLabels.erase(lookupFloat[0]->getLabel());
1074 //restore real lastlabel to save below
1075 lookupFloat[0]->setLabel(saveLabel);
1079 lastLabel = lookupFloat[0]->getLabel();
1081 //get next line to process
1082 //prevent memory leak
1083 for (int i = 0; i < lookupFloat.size(); i++) { delete lookupFloat[i]; }
1084 lookupFloat = input->getSharedRAbundFloatVectors();
1088 if (m->control_pressed) { delete input; return 0; }
1090 //output error messages about any remaining user labels
1091 set<string>::iterator it;
1092 bool needToRun = false;
1093 for (it = userLabels.begin(); it != userLabels.end(); it++) {
1094 m->mothurOut("Your file does not include the label " + *it);
1095 if (processedLabels.count(lastLabel) != 1) {
1096 m->mothurOut(". I will use " + lastLabel + "."); m->mothurOutEndLine();
1099 m->mothurOut(". Please refer to " + lastLabel + "."); m->mothurOutEndLine();
1103 //run last label if you need to
1104 if (needToRun == true) {
1105 for (int i = 0; i < lookupFloat.size(); i++) { if (lookupFloat[i] != NULL) { delete lookupFloat[i]; } }
1106 lookupFloat = input->getSharedRAbundFloatVectors(lastLabel);
1112 catch(exception& e) {
1113 m->errorOut(e, "IndicatorCommand", "getShared");
1117 //**********************************************************************************************************************
1118 vector<float> IndicatorCommand::driver(vector< vector<SharedRAbundFloatVector*> >& groupings, map< vector<int>, vector<int> > groupingsMap, int num, vector<float> indicatorValues, int numIters){
1120 vector<float> pvalues;
1121 pvalues.resize(indicatorValues.size(), 0);
1123 for(int i=0;i<numIters;i++){
1124 if (m->control_pressed) { break; }
1125 groupingsMap = randomizeGroupings(groupings, num);
1126 vector<float> randomIndicatorValues = getValues(groupings, groupingsMap);
1128 for (int j = 0; j < indicatorValues.size(); j++) {
1129 if (randomIndicatorValues[j] >= indicatorValues[j]) { pvalues[j]++; }
1135 }catch(exception& e) {
1136 m->errorOut(e, "IndicatorCommand", "driver");
1140 //**********************************************************************************************************************
1141 vector<float> IndicatorCommand::getPValues(vector< vector<SharedRAbundFloatVector*> >& groupings, map< vector<int>, vector<int> > groupingsMap, int num, vector<float> indicatorValues){
1143 vector<float> pvalues;
1145 #if defined (__APPLE__) || (__MACH__) || (linux) || (__linux) || (__linux__) || (__unix__) || (__unix)
1146 if(processors == 1){
1147 pvalues = driver(groupings, groupingsMap, num, indicatorValues, iters);
1148 for (int i = 0; i < pvalues.size(); i++) { pvalues[i] /= (double)iters; }
1151 //divide iters between processors
1152 int numItersPerProcessor = iters / processors;
1154 vector<int> processIDS;
1158 //loop through and create all the processes you want
1159 while (process != processors) {
1163 processIDS.push_back(pid); //create map from line number to pid so you can append files in correct order later
1165 }else if (pid == 0){
1166 pvalues = driver(groupings, groupingsMap, num, indicatorValues, numItersPerProcessor);
1168 //pass pvalues to parent
1170 string tempFile = toString(getpid()) + ".pvalues.temp";
1171 m->openOutputFile(tempFile, out);
1174 for (int i = 0; i < pvalues.size(); i++) {
1175 out << pvalues[i] << '\t';
1183 m->mothurOut("[ERROR]: unable to spawn the necessary processes."); m->mothurOutEndLine();
1184 for (int i = 0; i < processIDS.size(); i++) { kill (processIDS[i], SIGINT); }
1190 //special case for last processor in case it doesn't divide evenly
1191 numItersPerProcessor = iters - ((processors-1) * numItersPerProcessor);
1192 pvalues = driver(groupings, groupingsMap, num, indicatorValues, numItersPerProcessor);
1194 //force parent to wait until all the processes are done
1195 for (int i=0;i<processIDS.size();i++) {
1196 int temp = processIDS[i];
1201 for (int i = 0; i < processIDS.size(); i++) {
1203 string tempFile = toString(processIDS[i]) + ".pvalues.temp";
1204 m->openInputFile(tempFile, in);
1206 ////// to do ///////////
1207 int numTemp; numTemp = 0;
1208 for (int j = 0; j < pvalues.size(); j++) {
1209 in >> numTemp; m->gobble(in);
1210 pvalues[j] += numTemp;
1212 in.close(); m->mothurRemove(tempFile);
1214 for (int i = 0; i < pvalues.size(); i++) { pvalues[i] /= (double)iters; }
1217 pvalues = driver(groupings, groupingsMap, num, indicatorValues, iters);
1218 for (int i = 0; i < pvalues.size(); i++) { pvalues[i] /= (double)iters; }
1223 catch(exception& e) {
1224 m->errorOut(e, "IndicatorCommand", "getPValues");
1229 //**********************************************************************************************************************
1230 //same as above, just data type difference
1231 vector<float> IndicatorCommand::driver(vector< vector<SharedRAbundVector*> >& groupings, map< vector<int>, vector<int> > groupingsMap, int num, vector<float> indicatorValues, int numIters){
1233 vector<float> pvalues;
1234 pvalues.resize(indicatorValues.size(), 0);
1236 for(int i=0;i<numIters;i++){
1237 if (m->control_pressed) { break; }
1238 groupingsMap = randomizeGroupings(groupings, num);
1239 vector<float> randomIndicatorValues = getValues(groupings, groupingsMap);
1241 for (int j = 0; j < indicatorValues.size(); j++) {
1242 if (randomIndicatorValues[j] >= indicatorValues[j]) { pvalues[j]++; }
1248 }catch(exception& e) {
1249 m->errorOut(e, "IndicatorCommand", "driver");
1253 //**********************************************************************************************************************
1254 //same as above, just data type difference
1255 vector<float> IndicatorCommand::getPValues(vector< vector<SharedRAbundVector*> >& groupings, map< vector<int>, vector<int> > groupingsMap, int num, vector<float> indicatorValues){
1257 vector<float> pvalues;
1259 #if defined (__APPLE__) || (__MACH__) || (linux) || (__linux) || (__linux__) || (__unix__) || (__unix)
1260 if(processors == 1){
1261 pvalues = driver(groupings, groupingsMap, num, indicatorValues, iters);
1262 for (int i = 0; i < pvalues.size(); i++) { pvalues[i] /= (double)iters; }
1265 //divide iters between processors
1266 int numItersPerProcessor = iters / processors;
1268 vector<int> processIDS;
1271 //loop through and create all the processes you want
1272 while (process != processors) {
1276 processIDS.push_back(pid); //create map from line number to pid so you can append files in correct order later
1278 }else if (pid == 0){
1279 pvalues = driver(groupings, groupingsMap, num, indicatorValues, numItersPerProcessor);
1281 //pass pvalues to parent
1283 string tempFile = toString(getpid()) + ".pvalues.temp";
1284 m->openOutputFile(tempFile, out);
1287 for (int i = 0; i < pvalues.size(); i++) {
1288 out << pvalues[i] << '\t';
1296 m->mothurOut("[ERROR]: unable to spawn the necessary processes."); m->mothurOutEndLine();
1297 for (int i = 0; i < processIDS.size(); i++) { kill (processIDS[i], SIGINT); }
1303 //special case for last processor in case it doesn't divide evenly
1304 numItersPerProcessor = iters - ((processors-1) * numItersPerProcessor);
1305 pvalues = driver(groupings, groupingsMap, num, indicatorValues, numItersPerProcessor);
1307 //force parent to wait until all the processes are done
1308 for (int i=0;i<processIDS.size();i++) {
1309 int temp = processIDS[i];
1314 for (int i = 0; i < processIDS.size(); i++) {
1316 string tempFile = toString(processIDS[i]) + ".pvalues.temp";
1317 m->openInputFile(tempFile, in);
1319 ////// to do ///////////
1320 int numTemp; numTemp = 0;
1321 for (int j = 0; j < pvalues.size(); j++) {
1322 in >> numTemp; m->gobble(in);
1323 pvalues[j] += numTemp;
1325 in.close(); m->mothurRemove(tempFile);
1327 for (int i = 0; i < pvalues.size(); i++) { pvalues[i] /= (double)iters; }
1330 pvalues = driver(groupings, groupingsMap, num, indicatorValues, iters);
1331 for (int i = 0; i < pvalues.size(); i++) { pvalues[i] /= (double)iters; }
1336 catch(exception& e) {
1337 m->errorOut(e, "IndicatorCommand", "getPValues");
1341 //**********************************************************************************************************************
1342 //swap groups between groupings, in essence randomizing the second column of the design file
1343 map< vector<int>, vector<int> > IndicatorCommand::randomizeGroupings(vector< vector<SharedRAbundVector*> >& groupings, int numLookupGroups){
1346 map< vector<int>, vector<int> > randomGroupings;
1348 for (int i = 0; i < numLookupGroups; i++) {
1349 if (m->control_pressed) {break;}
1351 //get random groups to swap to switch with
1352 //generate random int between 0 and groupings.size()-1
1353 int z = m->getRandomIndex(groupings.size()-1);
1354 int x = m->getRandomIndex(groupings.size()-1);
1355 int a = m->getRandomIndex(groupings[z].size()-1);
1356 int b = m->getRandomIndex(groupings[x].size()-1);
1357 //cout << i << '\t' << z << '\t' << x << '\t' << a << '\t' << b << endl;
1358 //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; }
1363 from.push_back(z); from.push_back(a);
1364 to.push_back(x); to.push_back(b);
1366 randomGroupings[from] = to;
1368 //cout << "done" << endl;
1369 return randomGroupings;
1371 catch(exception& e) {
1372 m->errorOut(e, "IndicatorCommand", "randomizeGroupings");
1376 //**********************************************************************************************************************
1377 //swap groups between groupings, in essence randomizing the second column of the design file
1378 map< vector<int>, vector<int> > IndicatorCommand::randomizeGroupings(vector< vector<SharedRAbundFloatVector*> >& groupings, int numLookupGroups){
1381 map< vector<int>, vector<int> > randomGroupings;
1383 for (int i = 0; i < numLookupGroups; i++) {
1385 //get random groups to swap to switch with
1386 //generate random int between 0 and groupings.size()-1
1387 int z = m->getRandomIndex(groupings.size()-1);
1388 int x = m->getRandomIndex(groupings.size()-1);
1389 int a = m->getRandomIndex(groupings[z].size()-1);
1390 int b = m->getRandomIndex(groupings[x].size()-1);
1391 //cout << i << '\t' << z << '\t' << x << '\t' << a << '\t' << b << endl;
1396 from.push_back(z); from.push_back(a);
1397 to.push_back(x); to.push_back(b);
1399 randomGroupings[from] = to;
1402 return randomGroupings;
1404 catch(exception& e) {
1405 m->errorOut(e, "IndicatorCommand", "randomizeGroupings");
1410 /*****************************************************************/