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","summary",false,false,true); parameters.push_back(pdesign);
19 CommandParameter pshared("shared", "InputTypes", "", "", "SharedRel", "SharedRel", "none","summary",false,false,true); parameters.push_back(pshared);
20 CommandParameter prelabund("relabund", "InputTypes", "", "", "SharedRel", "SharedRel", "none","summary",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","tree-summary",false,false,true); 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::getOutputPattern(string type) {
65 if (type == "tree") { pattern = "[filename],indicator.tre"; }
66 else if (type == "summary") { pattern = "[filename],indicator.summary"; }
67 else { m->mothurOut("[ERROR]: No definition for type " + type + " output pattern.\n"); m->control_pressed = true; }
72 m->errorOut(e, "IndicatorCommand", "getOutputPattern");
76 //**********************************************************************************************************************
77 IndicatorCommand::IndicatorCommand(){
79 abort = true; calledHelp = true;
81 vector<string> tempOutNames;
82 outputTypes["tree"] = tempOutNames;
83 outputTypes["summary"] = tempOutNames;
86 m->errorOut(e, "IndicatorCommand", "IndicatorCommand");
90 //**********************************************************************************************************************
91 IndicatorCommand::IndicatorCommand(string option) {
93 abort = false; calledHelp = false;
95 //allow user to run help
96 if(option == "help") { help(); abort = true; calledHelp = true; }
97 else if(option == "citation") { citation(); abort = true; calledHelp = true;}
100 vector<string> myArray = setParameters();
102 OptionParser parser(option);
103 map<string, string> parameters = parser.getParameters();
105 ValidParameters validParameter;
106 map<string, string>::iterator it;
108 //check to make sure all parameters are valid for command
109 for (it = parameters.begin(); it != parameters.end(); it++) {
110 if (validParameter.isValidParameter(it->first, myArray, it->second) != true) { abort = true; }
116 m->Treenames.clear();
118 vector<string> tempOutNames;
119 outputTypes["tree"] = tempOutNames;
120 outputTypes["summary"] = tempOutNames;
122 //if the user changes the input directory command factory will send this info to us in the output parameter
123 string inputDir = validParameter.validFile(parameters, "inputdir", false);
124 if (inputDir == "not found"){ inputDir = ""; }
127 it = parameters.find("tree");
128 //user has given a template file
129 if(it != parameters.end()){
130 path = m->hasPath(it->second);
131 //if the user has not given a path then, add inputdir. else leave path alone.
132 if (path == "") { parameters["tree"] = inputDir + it->second; }
135 it = parameters.find("shared");
136 //user has given a template file
137 if(it != parameters.end()){
138 path = m->hasPath(it->second);
139 //if the user has not given a path then, add inputdir. else leave path alone.
140 if (path == "") { parameters["shared"] = inputDir + it->second; }
143 it = parameters.find("relabund");
144 //user has given a template file
145 if(it != parameters.end()){
146 path = m->hasPath(it->second);
147 //if the user has not given a path then, add inputdir. else leave path alone.
148 if (path == "") { parameters["relabund"] = inputDir + it->second; }
151 it = parameters.find("design");
152 //user has given a template file
153 if(it != parameters.end()){
154 path = m->hasPath(it->second);
155 //if the user has not given a path then, add inputdir. else leave path alone.
156 if (path == "") { parameters["design"] = inputDir + it->second; }
160 outputDir = validParameter.validFile(parameters, "outputdir", false); if (outputDir == "not found"){ outputDir = ""; }
162 //check for required parameters
163 treefile = validParameter.validFile(parameters, "tree", true);
164 if (treefile == "not open") { treefile = ""; abort = true; }
165 else if (treefile == "not found") { treefile = ""; }
166 else { m->setTreeFile(treefile); }
168 sharedfile = validParameter.validFile(parameters, "shared", true);
169 if (sharedfile == "not open") { abort = true; }
170 else if (sharedfile == "not found") { sharedfile = ""; }
171 else { inputFileName = sharedfile; m->setSharedFile(sharedfile); }
173 relabundfile = validParameter.validFile(parameters, "relabund", true);
174 if (relabundfile == "not open") { abort = true; }
175 else if (relabundfile == "not found") { relabundfile = ""; }
176 else { inputFileName = relabundfile; m->setRelAbundFile(relabundfile); }
178 designfile = validParameter.validFile(parameters, "design", true);
179 if (designfile == "not open") { designfile = ""; abort = true; }
180 else if (designfile == "not found") { designfile = ""; }
181 else { m->setDesignFile(designfile); }
183 groups = validParameter.validFile(parameters, "groups", false);
184 if (groups == "not found") { groups = ""; Groups.push_back("all"); }
185 else { m->splitAtDash(groups, Groups); }
186 m->setGroups(Groups);
188 label = validParameter.validFile(parameters, "label", false);
189 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=""; }
191 string temp = validParameter.validFile(parameters, "iters", false); if (temp == "not found") { temp = "1000"; }
192 m->mothurConvert(temp, iters);
194 temp = validParameter.validFile(parameters, "processors", false); if (temp == "not found"){ temp = m->getProcessors(); }
195 m->setProcessors(temp);
196 m->mothurConvert(temp, processors);
198 if ((relabundfile == "") && (sharedfile == "")) {
199 //is there are current file available for either of these?
200 //give priority to shared, then relabund
201 sharedfile = m->getSharedFile();
202 if (sharedfile != "") { inputFileName = sharedfile; m->mothurOut("Using " + sharedfile + " as input file for the shared parameter."); m->mothurOutEndLine(); }
204 relabundfile = m->getRelAbundFile();
205 if (relabundfile != "") { inputFileName = relabundfile; m->mothurOut("Using " + relabundfile + " as input file for the relabund parameter."); m->mothurOutEndLine(); }
207 m->mothurOut("No valid current files. You must provide a shared or relabund."); m->mothurOutEndLine();
214 if ((designfile == "") && (treefile == "")) {
215 treefile = m->getTreeFile();
216 if (treefile != "") { m->mothurOut("Using " + treefile + " as input file for the tree parameter."); m->mothurOutEndLine(); }
218 designfile = m->getDesignFile();
219 if (designfile != "") { m->mothurOut("Using " + designfile + " as input file for the design parameter."); m->mothurOutEndLine(); }
221 m->mothurOut("[ERROR]: You must provide either a tree or design file."); m->mothurOutEndLine(); abort = true;
226 if ((relabundfile != "") && (sharedfile != "")) { m->mothurOut("[ERROR]: You may not use both a shared and relabund file."); m->mothurOutEndLine(); abort = true; }
231 catch(exception& e) {
232 m->errorOut(e, "IndicatorCommand", "IndicatorCommand");
236 //**********************************************************************************************************************
238 int IndicatorCommand::execute(){
241 if (abort == true) { if (calledHelp) { return 0; } return 2; }
243 cout.setf(ios::fixed, ios::floatfield); cout.setf(ios::showpoint);
245 int start = time(NULL);
247 //read designfile if given and set up groups for read of sharedfiles
248 if (designfile != "") {
249 designMap = new GroupMap(designfile);
250 designMap->readDesignMap();
252 //fill Groups - checks for "all" and for any typo groups
254 vector<string> nameGroups = designMap->getNamesOfGroups();
255 util.setGroups(Groups, nameGroups);
256 designMap->setNamesOfGroups(nameGroups);
258 vector<string> namesSeqs = designMap->getNamesSeqs(Groups);
259 m->setGroups(namesSeqs);
262 /***************************************************/
263 // use smart distancing to get right sharedRabund //
264 /***************************************************/
265 if (sharedfile != "") {
267 if (m->control_pressed) { if (designfile != "") { delete designMap; } for (int i = 0; i < lookup.size(); i++) { delete lookup[i]; } return 0; }
268 if (lookup[0] == NULL) { m->mothurOut("[ERROR] reading shared file."); m->mothurOutEndLine(); return 0; }
271 if (m->control_pressed) { if (designfile != "") { delete designMap; } for (int i = 0; i < lookupFloat.size(); i++) { delete lookupFloat[i]; } return 0; }
272 if (lookupFloat[0] == NULL) { m->mothurOut("[ERROR] reading relabund file."); m->mothurOutEndLine(); return 0; }
275 //reset groups if needed
276 if (designfile != "") { m->setGroups(Groups); }
278 /***************************************************/
279 // reading tree info //
280 /***************************************************/
281 if (treefile != "") {
282 string groupfile = "";
283 m->setTreeFile(treefile);
284 Tree* tree = new Tree(treefile); delete tree; //extracts names from tree to make faked out groupmap
285 ct = new CountTable();
286 bool mismatch = false;
289 map<string, string> groupMap;
291 for (int i = 0; i < m->Treenames.size(); i++) {
292 nameMap.insert(m->Treenames[i]);
293 //sanity check - is this a group that is not in the sharedfile?
294 if (i == 0) { gps.insert("Group1"); }
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 groupMap[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 groupMap[m->Treenames[i]] = "Group1";
314 ct->createTable(nameMap, groupMap, gps);
316 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; }
318 if (mismatch) { //cleanup and exit
319 if (designfile != "") { delete designMap; }
320 if (sharedfile != "") { for (int i = 0; i < lookup.size(); i++) { delete lookup[i]; } }
321 else { for (int i = 0; i < lookupFloat.size(); i++) { delete lookupFloat[i]; } }
326 read = new ReadNewickTree(treefile);
327 int readOk = read->read(ct);
329 if (readOk != 0) { m->mothurOut("Read Terminated."); m->mothurOutEndLine(); delete ct; delete read; return 0; }
331 vector<Tree*> T = read->getTrees();
335 if (m->control_pressed) {
336 if (designfile != "") { delete designMap; }
337 if (sharedfile != "") { for (int i = 0; i < lookup.size(); i++) { delete lookup[i]; } }
338 else { for (int i = 0; i < lookupFloat.size(); i++) { delete lookupFloat[i]; } }
339 for (int i = 0; i < T.size(); i++) { delete T[i]; } delete ct; return 0;
342 T[0]->assembleTree();
344 /***************************************************/
345 // create ouptut tree - respecting pickedGroups //
346 /***************************************************/
347 Tree* outputTree = new Tree(m->getNumGroups(), ct);
349 outputTree->getSubTree(T[0], m->getGroups());
350 outputTree->assembleTree();
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 ct; return 0;
362 /***************************************************/
363 // get indicator species values //
364 /***************************************************/
365 GetIndicatorSpecies(outputTree);
366 delete outputTree; delete ct;
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 map<string, string> variables;
409 variables["[filename]"] = thisOutputDir + m->getRootName(m->getSimpleName(inputFileName));
410 string outputFileName = getOutputFileName("summary", variables);
411 outputNames.push_back(outputFileName); outputTypes["summary"].push_back(outputFileName);
414 m->openOutputFile(outputFileName, out);
415 out.setf(ios::fixed, ios::floatfield); out.setf(ios::showpoint);
416 m->mothurOutEndLine(); m->mothurOut("Species\tIndicator_Groups\tIndicatorValue\tpValue\n");
419 if (sharedfile != "") { numBins = lookup[0]->getNumBins(); }
420 else { numBins = lookupFloat[0]->getNumBins(); }
422 if (m->control_pressed) { out.close(); return 0; }
424 /*****************************************************/
425 //create vectors containing rabund info //
426 /*****************************************************/
428 vector<float> indicatorValues; //size of numBins
429 vector<float> pValues;
430 vector<string> indicatorGroups;
431 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.
433 if (sharedfile != "") {
434 vector< vector<SharedRAbundVector*> > groupings;
435 set<string> groupsAlreadyAdded;
436 vector<SharedRAbundVector*> subset;
439 for (int i = 0; i < (designMap->getNamesOfGroups()).size(); i++) {
441 for (int k = 0; k < lookup.size(); k++) {
442 //are you from this grouping?
443 if (designMap->getGroup(lookup[k]->getGroup()) == (designMap->getNamesOfGroups())[i]) {
444 subset.push_back(lookup[k]);
445 groupsAlreadyAdded.insert(lookup[k]->getGroup());
448 if (subset.size() != 0) { groupings.push_back(subset); }
452 if (groupsAlreadyAdded.size() != lookup.size()) { m->mothurOut("[ERROR]: could not make proper groupings."); m->mothurOutEndLine(); }
454 indicatorValues = getValues(groupings, indicatorGroups, randomGroupingsMap);
456 pValues = getPValues(groupings, randomGroupingsMap, lookup.size(), indicatorValues);
458 vector< vector<SharedRAbundFloatVector*> > groupings;
459 set<string> groupsAlreadyAdded;
460 vector<SharedRAbundFloatVector*> subset;
463 for (int i = 0; i < (designMap->getNamesOfGroups()).size(); i++) {
464 for (int k = 0; k < lookupFloat.size(); k++) {
465 //are you from this grouping?
466 if (designMap->getGroup(lookupFloat[k]->getGroup()) == (designMap->getNamesOfGroups())[i]) {
467 subset.push_back(lookupFloat[k]);
468 groupsAlreadyAdded.insert(lookupFloat[k]->getGroup());
471 if (subset.size() != 0) { groupings.push_back(subset); }
475 if (groupsAlreadyAdded.size() != lookupFloat.size()) { m->mothurOut("[ERROR]: could not make proper groupings."); m->mothurOutEndLine(); }
477 indicatorValues = getValues(groupings, indicatorGroups, randomGroupingsMap);
479 pValues = getPValues(groupings, randomGroupingsMap, lookupFloat.size(), indicatorValues);
482 if (m->control_pressed) { out.close(); return 0; }
485 /******************************************************/
486 //output indicator values to table form //
487 /*****************************************************/
488 out << "OTU\tIndicator_Groups\tIndicator_Value\tpValue" << endl;
489 for (int j = 0; j < indicatorValues.size(); j++) {
491 if (m->control_pressed) { out.close(); return 0; }
493 out << m->currentBinLabels[j] << '\t' << indicatorGroups[j] << '\t' << indicatorValues[j] << '\t';
495 if (pValues[j] > (1/(float)iters)) { out << pValues[j] << endl; }
496 else { out << "<" << (1/(float)iters) << endl; }
498 if (pValues[j] <= 0.05) {
499 cout << m->currentBinLabels[j] << '\t' << indicatorGroups[j] << '\t' << indicatorValues[j] << '\t';
500 string pValueString = "<" + toString((1/(float)iters));
501 if (pValues[j] > (1/(float)iters)) { pValueString = toString(pValues[j]); cout << pValues[j];}
502 else { cout << "<" << (1/(float)iters); }
503 m->mothurOutJustToLog(m->currentBinLabels[j] + "\t" + indicatorGroups[j] + "\t" + toString(indicatorValues[j]) + "\t" + pValueString);
504 m->mothurOutEndLine();
512 catch(exception& e) {
513 m->errorOut(e, "IndicatorCommand", "GetIndicatorSpecies");
517 //**********************************************************************************************************************
518 //traverse tree finding indicator species values for each otu at each node
519 //label node with otu number that has highest indicator value
520 //report all otu values to file
521 int IndicatorCommand::GetIndicatorSpecies(Tree*& T){
524 string thisOutputDir = outputDir;
525 if (outputDir == "") { thisOutputDir += m->hasPath(inputFileName); }
526 map<string, string> variables;
527 variables["[filename]"] = thisOutputDir + m->getRootName(m->getSimpleName(inputFileName));
528 string outputFileName = getOutputFileName("summary",variables);
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] << "_IndGroups" << '\t' << m->currentBinLabels[i] << "_IndValue" << '\t' << "pValue" << '\t'; }
544 m->mothurOutEndLine(); m->mothurOut("Node\tSpecies\tIndicator_Groups\tIndicatorValue\tpValue\n");
546 string treeOutputDir = outputDir;
547 if (outputDir == "") { treeOutputDir += m->hasPath(treefile); }
548 variables["[filename]"] = treeOutputDir + m->getRootName(m->getSimpleName(treefile));
549 string outputTreeFileName = getOutputFileName("tree", variables);
552 //create a map from tree node index to names of descendants, save time later to know which sharedRabund you need
553 map<int, set<string> > nodeToDescendants;
554 map<int, set<int> > descendantNodes;
555 for (int i = 0; i < T->getNumNodes(); i++) {
556 if (m->control_pressed) { return 0; }
558 nodeToDescendants[i] = getDescendantList(T, i, nodeToDescendants, descendantNodes);
561 //you need the distances to leaf to decide grouping below
562 //this will also set branch lengths if the tree does not include them
563 map<int, float> distToRoot = getDistToRoot(T);
566 for (int i = T->getNumLeaves(); i < T->getNumNodes(); i++) {
567 //cout << endl << i+1 << endl;
568 if (m->control_pressed) { out.close(); return 0; }
570 /*****************************************************/
571 //create vectors containing rabund info //
572 /*****************************************************/
574 vector<float> indicatorValues; //size of numBins
575 vector<float> pValues;
576 vector<string> indicatorGroups;
577 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.
579 if (sharedfile != "") {
580 vector< vector<SharedRAbundVector*> > groupings;
582 //get nodes that will be a valid grouping
583 //you are valid if you are not one of my descendants
584 //AND your distToRoot is >= mine
585 //AND you were not added as part of a larger grouping. Largest nodes are added first.
587 set<string> groupsAlreadyAdded;
588 //create a grouping with my grouping
589 vector<SharedRAbundVector*> subset;
591 int doneCount = nodeToDescendants[i].size();
592 for (int k = 0; k < lookup.size(); k++) {
593 //is this descendant of i
594 if ((nodeToDescendants[i].count(lookup[k]->getGroup()) != 0)) {
595 subset.push_back(lookup[k]);
596 groupsAlreadyAdded.insert(lookup[k]->getGroup());
599 if (count == doneCount) { break; } //quit once you get the rabunds for this grouping
601 if (subset.size() != 0) { groupings.push_back(subset); }
604 for (int j = (T->getNumNodes()-1); j >= 0; j--) {
607 if ((descendantNodes[i].count(j) == 0) && (distToRoot[j] >= distToRoot[i])) {
608 vector<SharedRAbundVector*> subset;
610 int doneCount = nodeToDescendants[j].size();
611 for (int k = 0; k < lookup.size(); k++) {
612 //is this descendant of j, and we didn't already add this as part of a larger grouping
613 if ((nodeToDescendants[j].count(lookup[k]->getGroup()) != 0) && (groupsAlreadyAdded.count(lookup[k]->getGroup()) == 0)) {
614 subset.push_back(lookup[k]);
615 groupsAlreadyAdded.insert(lookup[k]->getGroup());
618 if (count == doneCount) { break; } //quit once you get the rabunds for this grouping
621 //if subset.size == 0 then the node was added as part of a larger grouping
622 if (subset.size() != 0) { groupings.push_back(subset); }
626 if (groupsAlreadyAdded.size() != lookup.size()) { m->mothurOut("[ERROR]: could not make proper groupings."); m->mothurOutEndLine(); }
628 indicatorValues = getValues(groupings, indicatorGroups, randomGroupingsMap);
630 pValues = getPValues(groupings, randomGroupingsMap, lookup.size(), indicatorValues);
632 vector< vector<SharedRAbundFloatVector*> > groupings;
634 //get nodes that will be a valid grouping
635 //you are valid if you are not one of my descendants
636 //AND your distToRoot is >= mine
637 //AND you were not added as part of a larger grouping. Largest nodes are added first.
639 set<string> groupsAlreadyAdded;
640 //create a grouping with my grouping
641 vector<SharedRAbundFloatVector*> subset;
643 int doneCount = nodeToDescendants[i].size();
644 for (int k = 0; k < lookupFloat.size(); k++) {
645 //is this descendant of i
646 if ((nodeToDescendants[i].count(lookupFloat[k]->getGroup()) != 0)) {
647 subset.push_back(lookupFloat[k]);
648 groupsAlreadyAdded.insert(lookupFloat[k]->getGroup());
651 if (count == doneCount) { break; } //quit once you get the rabunds for this grouping
653 if (subset.size() != 0) { groupings.push_back(subset); }
655 for (int j = (T->getNumNodes()-1); j >= 0; j--) {
656 if ((descendantNodes[i].count(j) == 0) && (distToRoot[j] >= distToRoot[i])) {
657 vector<SharedRAbundFloatVector*> subset;
659 int doneCount = nodeToDescendants[j].size();
660 for (int k = 0; k < lookupFloat.size(); k++) {
661 //is this descendant of j, and we didn't already add this as part of a larger grouping
662 if ((nodeToDescendants[j].count(lookupFloat[k]->getGroup()) != 0) && (groupsAlreadyAdded.count(lookupFloat[k]->getGroup()) == 0)) {
663 subset.push_back(lookupFloat[k]);
664 groupsAlreadyAdded.insert(lookupFloat[k]->getGroup());
667 if (count == doneCount) { break; } //quit once you get the rabunds for this grouping
670 //if subset.size == 0 then the node was added as part of a larger grouping
671 if (subset.size() != 0) { groupings.push_back(subset); }
675 if (groupsAlreadyAdded.size() != lookupFloat.size()) { m->mothurOut("[ERROR]: could not make proper groupings."); m->mothurOutEndLine(); }
677 indicatorValues = getValues(groupings, indicatorGroups, randomGroupingsMap);
679 pValues = getPValues(groupings, randomGroupingsMap, lookupFloat.size(), indicatorValues);
682 if (m->control_pressed) { out.close(); return 0; }
685 /******************************************************/
686 //output indicator values to table form + label tree //
687 /*****************************************************/
688 out << (i+1) << '\t';
689 for (int j = 0; j < indicatorValues.size(); j++) {
691 if (m->control_pressed) { out.close(); return 0; }
693 if (pValues[j] < (1/(float)iters)) {
694 out << indicatorGroups[j] << '\t' << indicatorValues[j] << '\t' << '<' << (1/(float)iters) << '\t';
696 out << indicatorGroups[j] << '\t' << indicatorValues[j] << '\t' << pValues[j] << '\t';
699 if (pValues[j] <= 0.05) {
700 cout << i+1 << '\t' << m->currentBinLabels[j] << '\t' << indicatorGroups[j] << '\t' << indicatorValues[j] << '\t';
701 string pValueString = "<" + toString((1/(float)iters));
702 if (pValues[j] > (1/(float)iters)) { pValueString = toString(pValues[j]); cout << pValues[j];}
703 else { cout << "<" << (1/(float)iters); }
704 m->mothurOutJustToLog(toString(i) + "\t" + m->currentBinLabels[j] + "\t" + indicatorGroups[j] + "\t" + toString(indicatorValues[j]) + "\t" + pValueString);
705 m->mothurOutEndLine();
710 T->tree[i].setLabel((i+1));
715 m->openOutputFile(outputTreeFileName, outTree);
716 outputNames.push_back(outputTreeFileName); outputTypes["tree"].push_back(outputTreeFileName);
718 T->print(outTree, "both");
723 catch(exception& e) {
724 m->errorOut(e, "IndicatorCommand", "GetIndicatorSpecies");
728 //**********************************************************************************************************************
729 vector<float> IndicatorCommand::getValues(vector< vector<SharedRAbundFloatVector*> >& groupings, vector<string>& indicatorGroupings, map< vector<int>, vector<int> > groupingsMap){
731 vector<float> values;
732 map< vector<int>, vector<int> >::iterator it;
734 indicatorGroupings.clear();
736 //create grouping strings
737 vector<string> groupingsGroups;
738 for (int j = 0; j < groupings.size(); j++) {
739 string tempGrouping = "";
740 for (int k = 0; k < groupings[j].size()-1; k++) {
741 tempGrouping += groupings[j][k]->getGroup() + "-";
743 tempGrouping += groupings[j][groupings[j].size()-1]->getGroup();
744 groupingsGroups.push_back(tempGrouping);
749 for (int i = 0; i < groupings[0][0]->getNumBins(); i++) {
751 if (m->control_pressed) { return values; }
754 float AijDenominator = 0.0;
757 //get overall abundance of each grouping
758 for (int j = 0; j < groupings.size(); j++) {
760 float totalAbund = 0;
762 for (int k = 0; k < groupings[j].size(); k++) {
763 vector<int> temp; temp.push_back(j); temp.push_back(k);
764 it = groupingsMap.find(temp);
766 if (it == groupingsMap.end()) { //this one didnt get moved
767 totalAbund += groupings[j][k]->getAbundance(i);
768 if (groupings[j][k]->getAbundance(i) != 0.0) { numNotZero++; }
770 totalAbund += groupings[(it->second)[0]][(it->second)[1]]->getAbundance(i);
771 if (groupings[(it->second)[0]][(it->second)[1]]->getAbundance(i) != 0.0) { numNotZero++; }
777 float Aij = (totalAbund / (float) groupings[j].size());
778 terms.push_back(Aij);
780 //percentage of sites represented
781 Bij.push_back(numNotZero / (float) groupings[j].size());
783 AijDenominator += Aij;
786 float maxIndVal = 0.0;
787 string maxGrouping = "";
788 for (int j = 0; j < terms.size(); j++) {
789 float thisAij = (terms[j] / AijDenominator); //relative abundance
790 float thisValue = thisAij * Bij[j] * 100.0;
793 if (thisValue > maxIndVal) { maxIndVal = thisValue; maxGrouping = groupingsGroups[j]; }
796 values.push_back(maxIndVal);
797 indicatorGroupings.push_back(maxGrouping);
802 catch(exception& e) {
803 m->errorOut(e, "IndicatorCommand", "getValues");
807 //**********************************************************************************************************************
808 //same as above, just data type difference
809 vector<float> IndicatorCommand::getValues(vector< vector<SharedRAbundVector*> >& groupings, vector<string>& indicatorGroupings, map< vector<int>, vector<int> > groupingsMap){
811 vector<float> values;
812 map< vector<int>, vector<int> >::iterator it;
814 indicatorGroupings.clear();
816 //create grouping strings
817 vector<string> groupingsGroups;
818 for (int j = 0; j < groupings.size(); j++) {
819 string tempGrouping = "";
820 for (int k = 0; k < groupings[j].size()-1; k++) {
821 tempGrouping += groupings[j][k]->getGroup() + "-";
823 tempGrouping += groupings[j][groupings[j].size()-1]->getGroup();
824 groupingsGroups.push_back(tempGrouping);
829 for (int i = 0; i < groupings[0][0]->getNumBins(); i++) {
831 float AijDenominator = 0.0;
833 //get overall abundance of each grouping
834 for (int j = 0; j < groupings.size(); j++) {
836 int totalAbund = 0.0;
838 for (int k = 0; k < groupings[j].size(); k++) {
839 vector<int> temp; temp.push_back(j); temp.push_back(k);
840 it = groupingsMap.find(temp);
842 if (it == groupingsMap.end()) { //this one didnt get moved
843 totalAbund += groupings[j][k]->getAbundance(i);
844 if (groupings[j][k]->getAbundance(i) != 0.0) { numNotZero++; }
846 totalAbund += groupings[(it->second)[0]][(it->second)[1]]->getAbundance(i);
847 if (groupings[(it->second)[0]][(it->second)[1]]->getAbundance(i) != 0.0) { numNotZero++; }
853 float Aij = (totalAbund / (float) groupings[j].size());
854 terms.push_back(Aij);
856 //percentage of sites represented
857 Bij.push_back(numNotZero / (float) groupings[j].size());
859 AijDenominator += Aij;
862 float maxIndVal = 0.0;
863 string maxGrouping = "";
864 for (int j = 0; j < terms.size(); j++) {
865 float thisAij = (terms[j] / AijDenominator); //relative abundance
866 float thisValue = thisAij * Bij[j] * 100.0;
869 if (thisValue > maxIndVal) { maxIndVal = thisValue; maxGrouping = groupingsGroups[j]; }
872 values.push_back(maxIndVal);
873 indicatorGroupings.push_back(maxGrouping);
878 catch(exception& e) {
879 m->errorOut(e, "IndicatorCommand", "getValues");
883 //**********************************************************************************************************************
884 //you need the distances to root to decide groupings
885 //this will also set branch lengths if the tree does not include them
886 map<int, float> IndicatorCommand::getDistToRoot(Tree*& T){
888 map<int, float> dists;
890 bool hasBranchLengths = false;
891 for (int i = 0; i < T->getNumNodes(); i++) {
892 if (T->tree[i].getBranchLength() > 0.0) { hasBranchLengths = true; break; }
895 //set branchlengths if needed
896 if (!hasBranchLengths) {
897 for (int i = 0; i < T->getNumNodes(); i++) {
899 int lc = T->tree[i].getLChild();
900 int rc = T->tree[i].getRChild();
902 if (lc == -1) { // you are a leaf
903 //if you are a leaf set you priliminary length to 1.0, this may adjust later
904 T->tree[i].setBranchLength(1.0);
906 }else{ // you are an internal node
907 //look at your children's length to leaf
908 float ldist = dists[lc];
909 float rdist = dists[rc];
911 float greater = ldist;
912 if (rdist > greater) { greater = rdist; dists[i] = ldist + 1.0;}
913 else { dists[i] = rdist + 1.0; }
916 //branch length = difference + 1
917 T->tree[lc].setBranchLength((abs(ldist-greater) + 1.0));
918 T->tree[rc].setBranchLength((abs(rdist-greater) + 1.0));
925 for (int i = 0; i < T->getNumNodes(); i++) {
930 while(T->tree[index].getParent() != -1){
931 if (T->tree[index].getBranchLength() != -1) {
932 sum += abs(T->tree[index].getBranchLength());
934 index = T->tree[index].getParent();
942 catch(exception& e) {
943 m->errorOut(e, "IndicatorCommand", "getLengthToLeaf");
947 //**********************************************************************************************************************
948 set<string> IndicatorCommand::getDescendantList(Tree*& T, int i, map<int, set<string> > descendants, map<int, set<int> >& nodes){
952 set<string>::iterator it;
954 int lc = T->tree[i].getLChild();
955 int rc = T->tree[i].getRChild();
957 if (lc == -1) { //you are a leaf your only descendant is yourself
958 set<int> temp; temp.insert(i);
961 if (designfile == "") {
962 names.insert(T->tree[i].getName());
964 vector<string> myGroup; myGroup.push_back(T->tree[i].getName());
965 vector<string> myReps = designMap->getNamesSeqs(myGroup);
966 for (int k = 0; k < myReps.size(); k++) {
967 names.insert(myReps[k]);
971 }else{ //your descedants are the combination of your childrens descendants
972 names = descendants[lc];
973 nodes[i] = nodes[lc];
974 for (it = descendants[rc].begin(); it != descendants[rc].end(); it++) {
977 for (set<int>::iterator itNum = nodes[rc].begin(); itNum != nodes[rc].end(); itNum++) {
978 nodes[i].insert(*itNum);
980 //you are your own descendant
986 catch(exception& e) {
987 m->errorOut(e, "IndicatorCommand", "getDescendantList");
991 //**********************************************************************************************************************
992 int IndicatorCommand::getShared(){
994 InputData* input = new InputData(sharedfile, "sharedfile");
995 lookup = input->getSharedRAbundVectors();
996 string lastLabel = lookup[0]->getLabel();
998 if (label == "") { label = lastLabel; delete input; return 0; }
1000 //if the users enters label "0.06" and there is no "0.06" in their file use the next lowest label.
1001 set<string> labels; labels.insert(label);
1002 set<string> processedLabels;
1003 set<string> userLabels = labels;
1005 //as long as you are not at the end of the file or done wih the lines you want
1006 while((lookup[0] != NULL) && (userLabels.size() != 0)) {
1007 if (m->control_pressed) { delete input; return 0; }
1009 if(labels.count(lookup[0]->getLabel()) == 1){
1010 processedLabels.insert(lookup[0]->getLabel());
1011 userLabels.erase(lookup[0]->getLabel());
1015 if ((m->anyLabelsToProcess(lookup[0]->getLabel(), userLabels, "") == true) && (processedLabels.count(lastLabel) != 1)) {
1016 string saveLabel = lookup[0]->getLabel();
1018 for (int i = 0; i < lookup.size(); i++) { delete lookup[i]; }
1019 lookup = input->getSharedRAbundVectors(lastLabel);
1021 processedLabels.insert(lookup[0]->getLabel());
1022 userLabels.erase(lookup[0]->getLabel());
1024 //restore real lastlabel to save below
1025 lookup[0]->setLabel(saveLabel);
1029 lastLabel = lookup[0]->getLabel();
1031 //get next line to process
1032 //prevent memory leak
1033 for (int i = 0; i < lookup.size(); i++) { delete lookup[i]; }
1034 lookup = input->getSharedRAbundVectors();
1038 if (m->control_pressed) { delete input; return 0; }
1040 //output error messages about any remaining user labels
1041 set<string>::iterator it;
1042 bool needToRun = false;
1043 for (it = userLabels.begin(); it != userLabels.end(); it++) {
1044 m->mothurOut("Your file does not include the label " + *it);
1045 if (processedLabels.count(lastLabel) != 1) {
1046 m->mothurOut(". I will use " + lastLabel + "."); m->mothurOutEndLine();
1049 m->mothurOut(". Please refer to " + lastLabel + "."); m->mothurOutEndLine();
1053 //run last label if you need to
1054 if (needToRun == true) {
1055 for (int i = 0; i < lookup.size(); i++) { if (lookup[i] != NULL) { delete lookup[i]; } }
1056 lookup = input->getSharedRAbundVectors(lastLabel);
1062 catch(exception& e) {
1063 m->errorOut(e, "IndicatorCommand", "getShared");
1067 //**********************************************************************************************************************
1068 int IndicatorCommand::getSharedFloat(){
1070 InputData* input = new InputData(relabundfile, "relabund");
1071 lookupFloat = input->getSharedRAbundFloatVectors();
1072 string lastLabel = lookupFloat[0]->getLabel();
1074 if (label == "") { label = lastLabel; delete input; return 0; }
1076 //if the users enters label "0.06" and there is no "0.06" in their file use the next lowest label.
1077 set<string> labels; labels.insert(label);
1078 set<string> processedLabels;
1079 set<string> userLabels = labels;
1081 //as long as you are not at the end of the file or done wih the lines you want
1082 while((lookupFloat[0] != NULL) && (userLabels.size() != 0)) {
1084 if (m->control_pressed) { delete input; return 0; }
1086 if(labels.count(lookupFloat[0]->getLabel()) == 1){
1087 processedLabels.insert(lookupFloat[0]->getLabel());
1088 userLabels.erase(lookupFloat[0]->getLabel());
1092 if ((m->anyLabelsToProcess(lookupFloat[0]->getLabel(), userLabels, "") == true) && (processedLabels.count(lastLabel) != 1)) {
1093 string saveLabel = lookupFloat[0]->getLabel();
1095 for (int i = 0; i < lookupFloat.size(); i++) { delete lookupFloat[i]; }
1096 lookupFloat = input->getSharedRAbundFloatVectors(lastLabel);
1098 processedLabels.insert(lookupFloat[0]->getLabel());
1099 userLabels.erase(lookupFloat[0]->getLabel());
1101 //restore real lastlabel to save below
1102 lookupFloat[0]->setLabel(saveLabel);
1106 lastLabel = lookupFloat[0]->getLabel();
1108 //get next line to process
1109 //prevent memory leak
1110 for (int i = 0; i < lookupFloat.size(); i++) { delete lookupFloat[i]; }
1111 lookupFloat = input->getSharedRAbundFloatVectors();
1115 if (m->control_pressed) { delete input; return 0; }
1117 //output error messages about any remaining user labels
1118 set<string>::iterator it;
1119 bool needToRun = false;
1120 for (it = userLabels.begin(); it != userLabels.end(); it++) {
1121 m->mothurOut("Your file does not include the label " + *it);
1122 if (processedLabels.count(lastLabel) != 1) {
1123 m->mothurOut(". I will use " + lastLabel + "."); m->mothurOutEndLine();
1126 m->mothurOut(". Please refer to " + lastLabel + "."); m->mothurOutEndLine();
1130 //run last label if you need to
1131 if (needToRun == true) {
1132 for (int i = 0; i < lookupFloat.size(); i++) { if (lookupFloat[i] != NULL) { delete lookupFloat[i]; } }
1133 lookupFloat = input->getSharedRAbundFloatVectors(lastLabel);
1139 catch(exception& e) {
1140 m->errorOut(e, "IndicatorCommand", "getShared");
1144 //**********************************************************************************************************************
1145 vector<float> IndicatorCommand::driver(vector< vector<SharedRAbundFloatVector*> >& groupings, map< vector<int>, vector<int> > groupingsMap, int num, vector<float> indicatorValues, int numIters){
1147 vector<float> pvalues;
1148 pvalues.resize(indicatorValues.size(), 0);
1149 vector<string> notUsedGroupings; //we dont care about the grouping for the pvalues since they are randomized, but we need to pass the function something to make it work.
1151 for(int i=0;i<numIters;i++){
1152 if (m->control_pressed) { break; }
1153 groupingsMap = randomizeGroupings(groupings, num);
1154 vector<float> randomIndicatorValues = getValues(groupings, notUsedGroupings, groupingsMap);
1156 for (int j = 0; j < indicatorValues.size(); j++) {
1157 if (randomIndicatorValues[j] >= indicatorValues[j]) { pvalues[j]++; }
1163 }catch(exception& e) {
1164 m->errorOut(e, "IndicatorCommand", "driver");
1168 //**********************************************************************************************************************
1169 vector<float> IndicatorCommand::getPValues(vector< vector<SharedRAbundFloatVector*> >& groupings, map< vector<int>, vector<int> > groupingsMap, int num, vector<float> indicatorValues){
1171 vector<float> pvalues;
1173 #if defined (__APPLE__) || (__MACH__) || (linux) || (__linux) || (__linux__) || (__unix__) || (__unix)
1174 if(processors == 1){
1175 pvalues = driver(groupings, groupingsMap, num, indicatorValues, iters);
1176 for (int i = 0; i < pvalues.size(); i++) { pvalues[i] /= (double)iters; }
1179 //divide iters between processors
1180 int numItersPerProcessor = iters / processors;
1182 vector<int> processIDS;
1186 //loop through and create all the processes you want
1187 while (process != processors) {
1191 processIDS.push_back(pid); //create map from line number to pid so you can append files in correct order later
1193 }else if (pid == 0){
1194 pvalues = driver(groupings, groupingsMap, num, indicatorValues, numItersPerProcessor);
1196 //pass pvalues to parent
1198 string tempFile = toString(getpid()) + ".pvalues.temp";
1199 m->openOutputFile(tempFile, out);
1202 for (int i = 0; i < pvalues.size(); i++) {
1203 out << pvalues[i] << '\t';
1211 m->mothurOut("[ERROR]: unable to spawn the necessary processes."); m->mothurOutEndLine();
1212 for (int i = 0; i < processIDS.size(); i++) { kill (processIDS[i], SIGINT); }
1218 //special case for last processor in case it doesn't divide evenly
1219 numItersPerProcessor = iters - ((processors-1) * numItersPerProcessor);
1220 pvalues = driver(groupings, groupingsMap, num, indicatorValues, numItersPerProcessor);
1222 //force parent to wait until all the processes are done
1223 for (int i=0;i<processIDS.size();i++) {
1224 int temp = processIDS[i];
1229 for (int i = 0; i < processIDS.size(); i++) {
1231 string tempFile = toString(processIDS[i]) + ".pvalues.temp";
1232 m->openInputFile(tempFile, in);
1234 ////// to do ///////////
1235 int numTemp; numTemp = 0;
1236 for (int j = 0; j < pvalues.size(); j++) {
1237 in >> numTemp; m->gobble(in);
1238 pvalues[j] += numTemp;
1240 in.close(); m->mothurRemove(tempFile);
1242 for (int i = 0; i < pvalues.size(); i++) { pvalues[i] /= (double)iters; }
1245 pvalues = driver(groupings, groupingsMap, num, indicatorValues, iters);
1246 for (int i = 0; i < pvalues.size(); i++) { pvalues[i] /= (double)iters; }
1251 catch(exception& e) {
1252 m->errorOut(e, "IndicatorCommand", "getPValues");
1257 //**********************************************************************************************************************
1258 //same as above, just data type difference
1259 vector<float> IndicatorCommand::driver(vector< vector<SharedRAbundVector*> >& groupings, map< vector<int>, vector<int> > groupingsMap, int num, vector<float> indicatorValues, int numIters){
1261 vector<float> pvalues;
1262 pvalues.resize(indicatorValues.size(), 0);
1263 vector<string> notUsedGroupings; //we dont care about the grouping for the pvalues since they are randomized, but we need to pass the function something to make it work.
1265 for(int i=0;i<numIters;i++){
1266 if (m->control_pressed) { break; }
1267 groupingsMap = randomizeGroupings(groupings, num);
1268 vector<float> randomIndicatorValues = getValues(groupings, notUsedGroupings, groupingsMap);
1270 for (int j = 0; j < indicatorValues.size(); j++) {
1271 if (randomIndicatorValues[j] >= indicatorValues[j]) { pvalues[j]++; }
1277 }catch(exception& e) {
1278 m->errorOut(e, "IndicatorCommand", "driver");
1282 //**********************************************************************************************************************
1283 //same as above, just data type difference
1284 vector<float> IndicatorCommand::getPValues(vector< vector<SharedRAbundVector*> >& groupings, map< vector<int>, vector<int> > groupingsMap, int num, vector<float> indicatorValues){
1286 vector<float> pvalues;
1288 #if defined (__APPLE__) || (__MACH__) || (linux) || (__linux) || (__linux__) || (__unix__) || (__unix)
1289 if(processors == 1){
1290 pvalues = driver(groupings, groupingsMap, num, indicatorValues, iters);
1291 for (int i = 0; i < pvalues.size(); i++) { pvalues[i] /= (double)iters; }
1294 //divide iters between processors
1295 int numItersPerProcessor = iters / processors;
1297 vector<int> processIDS;
1300 //loop through and create all the processes you want
1301 while (process != processors) {
1305 processIDS.push_back(pid); //create map from line number to pid so you can append files in correct order later
1307 }else if (pid == 0){
1308 pvalues = driver(groupings, groupingsMap, num, indicatorValues, numItersPerProcessor);
1310 //pass pvalues to parent
1312 string tempFile = toString(getpid()) + ".pvalues.temp";
1313 m->openOutputFile(tempFile, out);
1316 for (int i = 0; i < pvalues.size(); i++) {
1317 out << pvalues[i] << '\t';
1325 m->mothurOut("[ERROR]: unable to spawn the necessary processes."); m->mothurOutEndLine();
1326 for (int i = 0; i < processIDS.size(); i++) { kill (processIDS[i], SIGINT); }
1332 //special case for last processor in case it doesn't divide evenly
1333 numItersPerProcessor = iters - ((processors-1) * numItersPerProcessor);
1334 pvalues = driver(groupings, groupingsMap, num, indicatorValues, numItersPerProcessor);
1336 //force parent to wait until all the processes are done
1337 for (int i=0;i<processIDS.size();i++) {
1338 int temp = processIDS[i];
1343 for (int i = 0; i < processIDS.size(); i++) {
1345 string tempFile = toString(processIDS[i]) + ".pvalues.temp";
1346 m->openInputFile(tempFile, in);
1348 ////// to do ///////////
1349 int numTemp; numTemp = 0;
1350 for (int j = 0; j < pvalues.size(); j++) {
1351 in >> numTemp; m->gobble(in);
1352 pvalues[j] += numTemp;
1354 in.close(); m->mothurRemove(tempFile);
1356 for (int i = 0; i < pvalues.size(); i++) { pvalues[i] /= (double)iters; }
1359 pvalues = driver(groupings, groupingsMap, num, indicatorValues, iters);
1360 for (int i = 0; i < pvalues.size(); i++) { pvalues[i] /= (double)iters; }
1365 catch(exception& e) {
1366 m->errorOut(e, "IndicatorCommand", "getPValues");
1370 //**********************************************************************************************************************
1371 //swap groups between groupings, in essence randomizing the second column of the design file
1372 map< vector<int>, vector<int> > IndicatorCommand::randomizeGroupings(vector< vector<SharedRAbundVector*> >& groupings, int numLookupGroups){
1375 map< vector<int>, vector<int> > randomGroupings;
1377 for (int i = 0; i < numLookupGroups; i++) {
1378 if (m->control_pressed) {break;}
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;
1387 //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; }
1392 from.push_back(z); from.push_back(a);
1393 to.push_back(x); to.push_back(b);
1395 randomGroupings[from] = to;
1397 //cout << "done" << endl;
1398 return randomGroupings;
1400 catch(exception& e) {
1401 m->errorOut(e, "IndicatorCommand", "randomizeGroupings");
1405 //**********************************************************************************************************************
1406 //swap groups between groupings, in essence randomizing the second column of the design file
1407 map< vector<int>, vector<int> > IndicatorCommand::randomizeGroupings(vector< vector<SharedRAbundFloatVector*> >& groupings, int numLookupGroups){
1410 map< vector<int>, vector<int> > randomGroupings;
1412 for (int i = 0; i < numLookupGroups; i++) {
1414 //get random groups to swap to switch with
1415 //generate random int between 0 and groupings.size()-1
1416 int z = m->getRandomIndex(groupings.size()-1);
1417 int x = m->getRandomIndex(groupings.size()-1);
1418 int a = m->getRandomIndex(groupings[z].size()-1);
1419 int b = m->getRandomIndex(groupings[x].size()-1);
1420 //cout << i << '\t' << z << '\t' << x << '\t' << a << '\t' << b << endl;
1425 from.push_back(z); from.push_back(a);
1426 to.push_back(x); to.push_back(b);
1428 randomGroupings[from] = to;
1431 return randomGroupings;
1433 catch(exception& e) {
1434 m->errorOut(e, "IndicatorCommand", "randomizeGroupings");
1439 /*****************************************************************/