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, 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, 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, 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, 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, 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 map< vector<int>, vector<int> > 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, int num, vector<float> indicatorValues){
1171 vector<float> pvalues;
1173 if(processors == 1){
1174 pvalues = driver(groupings, num, indicatorValues, iters);
1175 for (int i = 0; i < pvalues.size(); i++) { pvalues[i] /= (double)iters; }
1177 //divide iters between processors
1178 vector<int> procIters;
1179 int numItersPerProcessor = iters / processors;
1181 //divide iters between processes
1182 for (int h = 0; h < processors; h++) {
1183 if(h == processors - 1){ numItersPerProcessor = iters - h * numItersPerProcessor; }
1184 procIters.push_back(numItersPerProcessor);
1187 vector<int> processIDS;
1190 #if defined (__APPLE__) || (__MACH__) || (linux) || (__linux) || (__linux__) || (__unix__) || (__unix)
1192 //loop through and create all the processes you want
1193 while (process != processors) {
1197 processIDS.push_back(pid); //create map from line number to pid so you can append files in correct order later
1199 }else if (pid == 0){
1200 pvalues = driver(groupings, num, indicatorValues, procIters[process]);
1202 //pass pvalues to parent
1204 string tempFile = toString(getpid()) + ".pvalues.temp";
1205 m->openOutputFile(tempFile, out);
1208 for (int i = 0; i < pvalues.size(); i++) {
1209 out << pvalues[i] << '\t';
1217 m->mothurOut("[ERROR]: unable to spawn the necessary processes."); m->mothurOutEndLine();
1218 for (int i = 0; i < processIDS.size(); i++) { kill (processIDS[i], SIGINT); }
1224 pvalues = driver(groupings, num, indicatorValues, procIters[0]);
1226 //force parent to wait until all the processes are done
1227 for (int i=0;i<processIDS.size();i++) {
1228 int temp = processIDS[i];
1233 for (int i = 0; i < processIDS.size(); i++) {
1235 string tempFile = toString(processIDS[i]) + ".pvalues.temp";
1236 m->openInputFile(tempFile, in);
1238 ////// to do ///////////
1239 int numTemp; numTemp = 0;
1240 for (int j = 0; j < pvalues.size(); j++) {
1241 in >> numTemp; m->gobble(in);
1242 pvalues[j] += numTemp;
1244 in.close(); m->mothurRemove(tempFile);
1246 for (int i = 0; i < pvalues.size(); i++) { pvalues[i] /= (double)iters; }
1250 vector<indicatorData*> pDataArray;
1251 DWORD dwThreadIdArray[processors-1];
1252 HANDLE hThreadArray[processors-1];
1254 //Create processor worker threads.
1255 for( int i=1; i<processors; i++ ){
1257 //make copy of lookup so we don't get access violations
1258 vector< vector<SharedRAbundFloatVector*> > newGroupings;
1260 for (int k = 0; k < groupings.size(); k++) {
1261 vector<SharedRAbundFloatVector*> newLookup;
1262 for (int l = 0; l < groupings[k].size(); l++) {
1263 SharedRAbundFloatVector* temp = new SharedRAbundFloatVector();
1264 temp->setLabel(groupings[k][l]->getLabel());
1265 temp->setGroup(groupings[k][l]->getGroup());
1266 newLookup.push_back(temp);
1268 newGroupings.push_back(newLookup);
1272 for (int l = 0; l < groupings.size(); l++) {
1273 for (int k = 0; k < groupings[l][0]->getNumBins(); k++) {
1274 if (m->control_pressed) { for (int j = 0; j < newGroupings.size(); j++) { for (int u = 0; u < newGroupings[j].size(); u++) { delete newGroupings[j][u]; } } return pvalues; }
1276 for (int j = 0; j < groupings[l].size(); j++) { newGroupings[l][j]->push_back(groupings[l][j]->getAbundance(k), groupings[l][j]->getGroup()); }
1280 vector<float> copyIValues = indicatorValues;
1282 indicatorData* temp = new indicatorData(m, procIters[i], newGroupings, num, copyIValues);
1283 pDataArray.push_back(temp);
1284 processIDS.push_back(i);
1286 hThreadArray[i-1] = CreateThread(NULL, 0, MyIndicatorThreadFunction, pDataArray[i-1], 0, &dwThreadIdArray[i-1]);
1290 pvalues = driver(groupings, num, indicatorValues, procIters[0]);
1292 //Wait until all threads have terminated.
1293 WaitForMultipleObjects(processors-1, hThreadArray, TRUE, INFINITE);
1295 //Close all thread handles and free memory allocations.
1296 for(int i=0; i < pDataArray.size(); i++){
1297 for (int j = 0; j < pDataArray[i]->pvalues.size(); j++) { pvalues[j] += pDataArray[i]->pvalues[j]; }
1299 for (int l = 0; l < pDataArray[i]->groupings.size(); l++) {
1300 for (int j = 0; j < pDataArray[i]->groupings[l].size(); j++) { delete pDataArray[i]->groupings[l][j]; }
1303 CloseHandle(hThreadArray[i]);
1304 delete pDataArray[i];
1307 for (int i = 0; i < pvalues.size(); i++) { pvalues[i] /= (double)iters; }
1314 catch(exception& e) {
1315 m->errorOut(e, "IndicatorCommand", "getPValues");
1320 //**********************************************************************************************************************
1321 //same as above, just data type difference
1322 vector<float> IndicatorCommand::driver(vector< vector<SharedRAbundVector*> >& groupings, int num, vector<float> indicatorValues, int numIters){
1324 vector<float> pvalues;
1325 pvalues.resize(indicatorValues.size(), 0);
1326 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.
1328 for(int i=0;i<numIters;i++){
1329 if (m->control_pressed) { break; }
1330 map< vector<int>, vector<int> > groupingsMap = randomizeGroupings(groupings, num);
1331 vector<float> randomIndicatorValues = getValues(groupings, notUsedGroupings, groupingsMap);
1333 for (int j = 0; j < indicatorValues.size(); j++) {
1334 if (randomIndicatorValues[j] >= indicatorValues[j]) { pvalues[j]++; }
1340 }catch(exception& e) {
1341 m->errorOut(e, "IndicatorCommand", "driver");
1345 //**********************************************************************************************************************
1346 //same as above, just data type difference
1347 vector<float> IndicatorCommand::getPValues(vector< vector<SharedRAbundVector*> >& groupings, int num, vector<float> indicatorValues){
1349 vector<float> pvalues;
1351 if(processors == 1){
1352 pvalues = driver(groupings, num, indicatorValues, iters);
1353 for (int i = 0; i < pvalues.size(); i++) { pvalues[i] /= (double)iters; }
1355 //divide iters between processors
1356 vector<int> procIters;
1357 int numItersPerProcessor = iters / processors;
1359 //divide iters between processes
1360 for (int h = 0; h < processors; h++) {
1361 if(h == processors - 1){ numItersPerProcessor = iters - h * numItersPerProcessor; }
1362 procIters.push_back(numItersPerProcessor);
1365 vector<int> processIDS;
1368 #if defined (__APPLE__) || (__MACH__) || (linux) || (__linux) || (__linux__) || (__unix__) || (__unix)
1370 //loop through and create all the processes you want
1371 while (process != processors) {
1375 processIDS.push_back(pid); //create map from line number to pid so you can append files in correct order later
1377 }else if (pid == 0){
1378 pvalues = driver(groupings, num, indicatorValues, procIters[process]);
1380 //pass pvalues to parent
1382 string tempFile = toString(getpid()) + ".pvalues.temp";
1383 m->openOutputFile(tempFile, out);
1386 for (int i = 0; i < pvalues.size(); i++) {
1387 out << pvalues[i] << '\t';
1395 m->mothurOut("[ERROR]: unable to spawn the necessary processes."); m->mothurOutEndLine();
1396 for (int i = 0; i < processIDS.size(); i++) { kill (processIDS[i], SIGINT); }
1402 pvalues = driver(groupings, num, indicatorValues, procIters[0]);
1404 //force parent to wait until all the processes are done
1405 for (int i=0;i<processIDS.size();i++) {
1406 int temp = processIDS[i];
1411 for (int i = 0; i < processIDS.size(); i++) {
1413 string tempFile = toString(processIDS[i]) + ".pvalues.temp";
1414 m->openInputFile(tempFile, in);
1416 ////// to do ///////////
1417 int numTemp; numTemp = 0;
1418 for (int j = 0; j < pvalues.size(); j++) {
1419 in >> numTemp; m->gobble(in);
1420 pvalues[j] += numTemp;
1422 in.close(); m->mothurRemove(tempFile);
1424 for (int i = 0; i < pvalues.size(); i++) { pvalues[i] /= (double)iters; }
1428 vector<indicatorData*> pDataArray;
1429 DWORD dwThreadIdArray[processors-1];
1430 HANDLE hThreadArray[processors-1];
1432 //Create processor worker threads.
1433 for( int i=1; i<processors; i++ ){
1435 //make copy of lookup so we don't get access violations
1436 vector< vector<SharedRAbundFloatVector*> > newGroupings;
1438 for (int k = 0; k < groupings.size(); k++) {
1439 vector<SharedRAbundFloatVector*> newLookup;
1440 for (int l = 0; l < groupings[k].size(); l++) {
1441 SharedRAbundFloatVector* temp = new SharedRAbundFloatVector();
1442 temp->setLabel(groupings[k][l]->getLabel());
1443 temp->setGroup(groupings[k][l]->getGroup());
1444 newLookup.push_back(temp);
1446 newGroupings.push_back(newLookup);
1450 for (int l = 0; l < groupings.size(); l++) {
1451 for (int k = 0; k < groupings[l][0]->getNumBins(); k++) {
1452 if (m->control_pressed) { for (int j = 0; j < newGroupings.size(); j++) { for (int u = 0; u < newGroupings[j].size(); u++) { delete newGroupings[j][u]; } } return pvalues; }
1454 for (int j = 0; j < groupings[l].size(); j++) { newGroupings[l][j]->push_back((float)(groupings[l][j]->getAbundance(k)), groupings[l][j]->getGroup()); }
1458 vector<float> copyIValues = indicatorValues;
1460 indicatorData* temp = new indicatorData(m, procIters[i], newGroupings, num, copyIValues);
1461 pDataArray.push_back(temp);
1462 processIDS.push_back(i);
1464 hThreadArray[i-1] = CreateThread(NULL, 0, MyIndicatorThreadFunction, pDataArray[i-1], 0, &dwThreadIdArray[i-1]);
1468 pvalues = driver(groupings, num, indicatorValues, procIters[0]);
1470 //Wait until all threads have terminated.
1471 WaitForMultipleObjects(processors-1, hThreadArray, TRUE, INFINITE);
1473 //Close all thread handles and free memory allocations.
1474 for(int i=0; i < pDataArray.size(); i++){
1475 for (int j = 0; j < pDataArray[i]->pvalues.size(); j++) { pvalues[j] += pDataArray[i]->pvalues[j]; }
1477 for (int l = 0; l < pDataArray[i]->groupings.size(); l++) {
1478 for (int j = 0; j < pDataArray[i]->groupings[l].size(); j++) { delete pDataArray[i]->groupings[l][j]; }
1481 CloseHandle(hThreadArray[i]);
1482 delete pDataArray[i];
1485 for (int i = 0; i < pvalues.size(); i++) { pvalues[i] /= (double)iters; }
1492 catch(exception& e) {
1493 m->errorOut(e, "IndicatorCommand", "getPValues");
1497 //**********************************************************************************************************************
1498 //swap groups between groupings, in essence randomizing the second column of the design file
1499 map< vector<int>, vector<int> > IndicatorCommand::randomizeGroupings(vector< vector<SharedRAbundVector*> >& groupings, int numLookupGroups){
1502 map< vector<int>, vector<int> > randomGroupings;
1504 for (int i = 0; i < numLookupGroups; i++) {
1505 if (m->control_pressed) {break;}
1507 //get random groups to swap to switch with
1508 //generate random int between 0 and groupings.size()-1
1509 int z = m->getRandomIndex(groupings.size()-1);
1510 int x = m->getRandomIndex(groupings.size()-1);
1511 int a = m->getRandomIndex(groupings[z].size()-1);
1512 int b = m->getRandomIndex(groupings[x].size()-1);
1513 //cout << i << '\t' << z << '\t' << x << '\t' << a << '\t' << b << endl;
1514 //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; }
1519 from.push_back(z); from.push_back(a);
1520 to.push_back(x); to.push_back(b);
1522 randomGroupings[from] = to;
1524 //cout << "done" << endl;
1525 return randomGroupings;
1527 catch(exception& e) {
1528 m->errorOut(e, "IndicatorCommand", "randomizeGroupings");
1532 //**********************************************************************************************************************
1533 //swap groups between groupings, in essence randomizing the second column of the design file
1534 map< vector<int>, vector<int> > IndicatorCommand::randomizeGroupings(vector< vector<SharedRAbundFloatVector*> >& groupings, int numLookupGroups){
1537 map< vector<int>, vector<int> > randomGroupings;
1539 for (int i = 0; i < numLookupGroups; i++) {
1541 //get random groups to swap to switch with
1542 //generate random int between 0 and groupings.size()-1
1543 int z = m->getRandomIndex(groupings.size()-1);
1544 int x = m->getRandomIndex(groupings.size()-1);
1545 int a = m->getRandomIndex(groupings[z].size()-1);
1546 int b = m->getRandomIndex(groupings[x].size()-1);
1547 //cout << i << '\t' << z << '\t' << x << '\t' << a << '\t' << b << endl;
1552 from.push_back(z); from.push_back(a);
1553 to.push_back(x); to.push_back(b);
1555 randomGroupings[from] = to;
1558 return randomGroupings;
1560 catch(exception& e) {
1561 m->errorOut(e, "IndicatorCommand", "randomizeGroupings");
1566 /*****************************************************************/