X-Git-Url: https://git.donarmstrong.com/?p=mothur.git;a=blobdiff_plain;f=indicatorcommand.cpp;h=fd818acb44c4c5280c6a54971d3d39413307861d;hp=01627c0ec1b9aa5d7971a564f131924b179ffc94;hb=050a3ff02473a3d4c0980964e1a9ebe52e55d6b8;hpb=3247d888e7aafc4a65ec9062a94dfd166c2c5b1d diff --git a/indicatorcommand.cpp b/indicatorcommand.cpp index 01627c0..fd818ac 100644 --- a/indicatorcommand.cpp +++ b/indicatorcommand.cpp @@ -8,35 +8,76 @@ */ #include "indicatorcommand.h" +#include "sharedutilities.h" + + //********************************************************************************************************************** -vector IndicatorCommand::getValidParameters(){ +vector IndicatorCommand::setParameters(){ try { - string Array[] = {"tree","shared","relabund","label","groups","outputdir","inputdir"}; - vector myArray (Array, Array+(sizeof(Array)/sizeof(string))); + CommandParameter piters("iters", "Number", "", "1000", "", "", "","",false,false); parameters.push_back(piters); + CommandParameter pdesign("design", "InputTypes", "", "", "TreeDesign", "TreeDesign", "none","summary",false,false,true); parameters.push_back(pdesign); + CommandParameter pshared("shared", "InputTypes", "", "", "SharedRel", "SharedRel", "none","summary",false,false,true); parameters.push_back(pshared); + CommandParameter prelabund("relabund", "InputTypes", "", "", "SharedRel", "SharedRel", "none","summary",false,false); parameters.push_back(prelabund); + CommandParameter pgroups("groups", "String", "", "", "", "", "","",false,false); parameters.push_back(pgroups); + CommandParameter plabel("label", "String", "", "", "", "", "","",false,false); parameters.push_back(plabel); + CommandParameter ptree("tree", "InputTypes", "", "", "TreeDesign", "TreeDesign", "none","tree-summary",false,false,true); parameters.push_back(ptree); + CommandParameter pinputdir("inputdir", "String", "", "", "", "", "","",false,false); parameters.push_back(pinputdir); + CommandParameter poutputdir("outputdir", "String", "", "", "", "", "","",false,false); parameters.push_back(poutputdir); + CommandParameter pprocessors("processors", "Number", "", "1", "", "", "","",false,false); parameters.push_back(pprocessors); + + vector myArray; + for (int i = 0; i < parameters.size(); i++) { myArray.push_back(parameters[i].name); } return myArray; } catch(exception& e) { - m->errorOut(e, "IndicatorCommand", "getValidParameters"); + m->errorOut(e, "IndicatorCommand", "setParameters"); exit(1); } } //********************************************************************************************************************** -vector IndicatorCommand::getRequiredParameters(){ +string IndicatorCommand::getHelpString(){ try { - string Array[] = {"tree"}; - vector myArray (Array, Array+(sizeof(Array)/sizeof(string))); - return myArray; + string helpString = ""; + 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"; + helpString += "The indicator command outputs a .indicator.summary file and a .indicator.tre if a tree is given. \n"; + 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"; + helpString += "The summary file lists the indicator value for each OTU for each node.\n"; + helpString += "The indicator command parameters are tree, groups, shared, relabund, design and label. \n"; + 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"; + 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"; + helpString += "The label parameter indicates at what distance your tree relates to the shared or relabund.\n"; + helpString += "The processors parameter allows you to specify how many processors you would like to use. The default is 1. \n"; + helpString += "The iters parameter allows you to set number of randomization for the P value. The default is 1000."; + helpString += "The indicator command should be used in the following format: indicator(tree=test.tre, shared=test.shared, label=0.03)\n"; + helpString += "Note: No spaces between parameter labels (i.e. tree), '=' and parameters (i.e.yourTreefile).\n"; + return helpString; } catch(exception& e) { - m->errorOut(e, "IndicatorCommand", "getRequiredParameters"); + m->errorOut(e, "IndicatorCommand", "getHelpString"); exit(1); } } //********************************************************************************************************************** +string IndicatorCommand::getOutputPattern(string type) { + try { + string pattern = ""; + + if (type == "tree") { pattern = "[filename],indicator.tre"; } + else if (type == "summary") { pattern = "[filename],indicator.summary"; } + else { m->mothurOut("[ERROR]: No definition for type " + type + " output pattern.\n"); m->control_pressed = true; } + + return pattern; + } + catch(exception& e) { + m->errorOut(e, "IndicatorCommand", "getOutputPattern"); + exit(1); + } +} +//********************************************************************************************************************** IndicatorCommand::IndicatorCommand(){ try { - abort = true; - //initialize outputTypes + abort = true; calledHelp = true; + setParameters(); vector tempOutNames; outputTypes["tree"] = tempOutNames; outputTypes["summary"] = tempOutNames; @@ -46,31 +87,17 @@ IndicatorCommand::IndicatorCommand(){ exit(1); } } - -//********************************************************************************************************************** -vector IndicatorCommand::getRequiredFiles(){ - try { - vector myArray; - return myArray; - } - catch(exception& e) { - m->errorOut(e, "IndicatorCommand", "getRequiredFiles"); - exit(1); - } -} //********************************************************************************************************************** IndicatorCommand::IndicatorCommand(string option) { try { - globaldata = GlobalData::getInstance(); - abort = false; + abort = false; calledHelp = false; //allow user to run help - if(option == "help") { help(); abort = true; } + if(option == "help") { help(); abort = true; calledHelp = true; } + else if(option == "citation") { citation(); abort = true; calledHelp = true;} else { - //valid paramters for this command - string Array[] = {"tree","shared","relabund","groups","label","outputdir","inputdir"}; - vector myArray (Array, Array+(sizeof(Array)/sizeof(string))); + vector myArray = setParameters(); OptionParser parser(option); map parameters = parser.getParameters(); @@ -83,7 +110,10 @@ IndicatorCommand::IndicatorCommand(string option) { if (validParameter.isValidParameter(it->first, myArray, it->second) != true) { abort = true; } } - globaldata->newRead(); + m->runParse = true; + m->clearGroups(); + m->clearAllGroups(); + m->Treenames.clear(); vector tempOutNames; outputTypes["tree"] = tempOutNames; @@ -118,37 +148,83 @@ IndicatorCommand::IndicatorCommand(string option) { if (path == "") { parameters["relabund"] = inputDir + it->second; } } + it = parameters.find("design"); + //user has given a template file + if(it != parameters.end()){ + path = m->hasPath(it->second); + //if the user has not given a path then, add inputdir. else leave path alone. + if (path == "") { parameters["design"] = inputDir + it->second; } + } } outputDir = validParameter.validFile(parameters, "outputdir", false); if (outputDir == "not found"){ outputDir = ""; } //check for required parameters treefile = validParameter.validFile(parameters, "tree", true); - if (treefile == "not open") { abort = true; } - else if (treefile == "not found") { treefile = ""; m->mothurOut("tree is a required parameter for the indicator command."); m->mothurOutEndLine(); abort = true; } - else { globaldata->setTreeFile(treefile); globaldata->setFormat("tree"); } + if (treefile == "not open") { treefile = ""; abort = true; } + else if (treefile == "not found") { treefile = ""; } + else { m->setTreeFile(treefile); } sharedfile = validParameter.validFile(parameters, "shared", true); if (sharedfile == "not open") { abort = true; } else if (sharedfile == "not found") { sharedfile = ""; } - else { inputFileName = sharedfile; } + else { inputFileName = sharedfile; m->setSharedFile(sharedfile); } relabundfile = validParameter.validFile(parameters, "relabund", true); if (relabundfile == "not open") { abort = true; } else if (relabundfile == "not found") { relabundfile = ""; } - else { inputFileName = relabundfile; } + else { inputFileName = relabundfile; m->setRelAbundFile(relabundfile); } + + designfile = validParameter.validFile(parameters, "design", true); + if (designfile == "not open") { designfile = ""; abort = true; } + else if (designfile == "not found") { designfile = ""; } + else { m->setDesignFile(designfile); } groups = validParameter.validFile(parameters, "groups", false); if (groups == "not found") { groups = ""; Groups.push_back("all"); } else { m->splitAtDash(groups, Groups); } - globaldata->Groups = Groups; + m->setGroups(Groups); label = validParameter.validFile(parameters, "label", false); 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=""; } - if ((relabundfile == "") && (sharedfile == "")) { m->mothurOut("You must provide either a shared or relabund file."); m->mothurOutEndLine(); abort = true; } + string temp = validParameter.validFile(parameters, "iters", false); if (temp == "not found") { temp = "1000"; } + m->mothurConvert(temp, iters); + + temp = validParameter.validFile(parameters, "processors", false); if (temp == "not found"){ temp = m->getProcessors(); } + m->setProcessors(temp); + m->mothurConvert(temp, processors); + + if ((relabundfile == "") && (sharedfile == "")) { + //is there are current file available for either of these? + //give priority to shared, then relabund + sharedfile = m->getSharedFile(); + if (sharedfile != "") { inputFileName = sharedfile; m->mothurOut("Using " + sharedfile + " as input file for the shared parameter."); m->mothurOutEndLine(); } + else { + relabundfile = m->getRelAbundFile(); + if (relabundfile != "") { inputFileName = relabundfile; m->mothurOut("Using " + relabundfile + " as input file for the relabund parameter."); m->mothurOutEndLine(); } + else { + m->mothurOut("No valid current files. You must provide a shared or relabund."); m->mothurOutEndLine(); + abort = true; + } + } + } + + + if ((designfile == "") && (treefile == "")) { + treefile = m->getTreeFile(); + if (treefile != "") { m->mothurOut("Using " + treefile + " as input file for the tree parameter."); m->mothurOutEndLine(); } + else { + designfile = m->getDesignFile(); + if (designfile != "") { m->mothurOut("Using " + designfile + " as input file for the design parameter."); m->mothurOutEndLine(); } + else { + m->mothurOut("[ERROR]: You must provide either a tree or design file."); m->mothurOutEndLine(); abort = true; + } + } + } + + if ((relabundfile != "") && (sharedfile != "")) { m->mothurOut("[ERROR]: You may not use both a shared and relabund file."); m->mothurOutEndLine(); abort = true; } - if ((relabundfile != "") && (sharedfile != "")) { m->mothurOut("You may not use both a shared and relabund file."); m->mothurOutEndLine(); abort = true; } } } @@ -159,121 +235,157 @@ IndicatorCommand::IndicatorCommand(string option) { } //********************************************************************************************************************** -void IndicatorCommand::help(){ - try { - m->mothurOut("The indicator command reads a shared or relabund file and a tree file, and outputs a .indicator.tre and .indicator.summary file. \n"); - m->mothurOut("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"); - m->mothurOut("The summary file lists the indicator value for each OTU for each node.\n"); - m->mothurOut("The indicator command parameters are tree, groups, shared, relabund and label. The tree parameter is required as well as either shared or relabund.\n"); - m->mothurOut("The groups parameter allows you to specify which of the groups in your shared or relabund you would like analyzed. The groups may be entered separated by dashes.\n"); - m->mothurOut("The label parameter indicates at what distance your tree relates to the shared or relabund.\n"); - m->mothurOut("The indicator command should be used in the following format: indicator(tree=test.tre, shared=test.shared, label=0.03)\n"); - m->mothurOut("Note: No spaces between parameter labels (i.e. tree), '=' and parameters (i.e.yourTreefile).\n\n"); - } - catch(exception& e) { - m->errorOut(e, "IndicatorCommand", "help"); - exit(1); - } -} - -//********************************************************************************************************************** - -IndicatorCommand::~IndicatorCommand(){} - -//********************************************************************************************************************** - int IndicatorCommand::execute(){ try { - if (abort == true) { return 0; } + if (abort == true) { if (calledHelp) { return 0; } return 2; } + + cout.setf(ios::fixed, ios::floatfield); cout.setf(ios::showpoint); + + int start = time(NULL); + + //read designfile if given and set up groups for read of sharedfiles + if (designfile != "") { + designMap = new GroupMap(designfile); + designMap->readDesignMap(); + + //fill Groups - checks for "all" and for any typo groups + SharedUtil util; + vector nameGroups = designMap->getNamesOfGroups(); + util.setGroups(Groups, nameGroups); + designMap->setNamesOfGroups(nameGroups); + + vector namesSeqs = designMap->getNamesSeqs(Groups); + m->setGroups(namesSeqs); + } /***************************************************/ // use smart distancing to get right sharedRabund // /***************************************************/ if (sharedfile != "") { getShared(); - if (m->control_pressed) { for (int i = 0; i < lookup.size(); i++) { delete lookup[i]; } return 0; } + if (m->control_pressed) { if (designfile != "") { delete designMap; } for (int i = 0; i < lookup.size(); i++) { delete lookup[i]; } return 0; } if (lookup[0] == NULL) { m->mothurOut("[ERROR] reading shared file."); m->mothurOutEndLine(); return 0; } }else { getSharedFloat(); - if (m->control_pressed) { for (int i = 0; i < lookupFloat.size(); i++) { delete lookupFloat[i]; } return 0; } + if (m->control_pressed) { if (designfile != "") { delete designMap; } for (int i = 0; i < lookupFloat.size(); i++) { delete lookupFloat[i]; } return 0; } if (lookupFloat[0] == NULL) { m->mothurOut("[ERROR] reading relabund file."); m->mothurOutEndLine(); return 0; } } + + //reset groups if needed + if (designfile != "") { m->setGroups(Groups); } /***************************************************/ // reading tree info // /***************************************************/ - string groupfile = ""; - Tree* tree = new Tree(treefile); delete tree; //extracts names from tree to make faked out groupmap - - globaldata->setGroupFile(groupfile); - treeMap = new TreeMap(); - bool mismatch = false; - for (int i = 0; i < globaldata->Treenames.size(); i++) { - //sanity check - is this a group that is not in the sharedfile? - if (!(m->inUsersGroups(globaldata->Treenames[i], globaldata->gGroupmap->namesOfGroups))) { - m->mothurOut("[ERROR]: " + globaldata->Treenames[i] + " is not a group in your shared or relabund file."); m->mothurOutEndLine(); - mismatch = true; - } - treeMap->addSeq(globaldata->Treenames[i], "Group1"); - } - - if (mismatch) { //cleanup and exit - if (sharedfile != "") { for (int i = 0; i < lookup.size(); i++) { delete lookup[i]; } } - else { for (int i = 0; i < lookupFloat.size(); i++) { delete lookupFloat[i]; } } - delete treeMap; - return 0; - } + if (treefile != "") { + string groupfile = ""; + m->setTreeFile(treefile); + Tree* tree = new Tree(treefile); delete tree; //extracts names from tree to make faked out groupmap + ct = new CountTable(); + bool mismatch = false; + + set nameMap; + map groupMap; + set gps; + for (int i = 0; i < m->Treenames.size(); i++) { + nameMap.insert(m->Treenames[i]); + //sanity check - is this a group that is not in the sharedfile? + if (i == 0) { gps.insert("Group1"); } + if (designfile == "") { + if (!(m->inUsersGroups(m->Treenames[i], m->getAllGroups()))) { + m->mothurOut("[ERROR]: " + m->Treenames[i] + " is not a group in your shared or relabund file."); m->mothurOutEndLine(); + mismatch = true; + } + groupMap[m->Treenames[i]] = "Group1"; + }else{ + vector myGroups; myGroups.push_back(m->Treenames[i]); + vector myNames = designMap->getNamesSeqs(myGroups); + + for(int k = 0; k < myNames.size(); k++) { + if (!(m->inUsersGroups(myNames[k], m->getAllGroups()))) { + m->mothurOut("[ERROR]: " + myNames[k] + " is not a group in your shared or relabund file."); m->mothurOutEndLine(); + mismatch = true; + } + } + groupMap[m->Treenames[i]] = "Group1"; + } + } + ct->createTable(nameMap, groupMap, gps); - globaldata->gTreemap = treeMap; - - read = new ReadNewickTree(treefile); - int readOk = read->read(); - - if (readOk != 0) { m->mothurOut("Read Terminated."); m->mothurOutEndLine(); globaldata->gTree.clear(); delete globaldata->gTreemap; delete read; return 0; } - - vector T = globaldata->gTree; - - delete read; - - if (m->control_pressed) { - if (sharedfile != "") { for (int i = 0; i < lookup.size(); i++) { delete lookup[i]; } } - else { for (int i = 0; i < lookupFloat.size(); i++) { delete lookupFloat[i]; } } - for (int i = 0; i < T.size(); i++) { delete T[i]; } globaldata->gTree.clear(); delete globaldata->gTreemap; return 0; - } + 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; } + + if (mismatch) { //cleanup and exit + if (designfile != "") { delete designMap; } + if (sharedfile != "") { for (int i = 0; i < lookup.size(); i++) { delete lookup[i]; } } + else { for (int i = 0; i < lookupFloat.size(); i++) { delete lookupFloat[i]; } } + delete ct; + return 0; + } + + read = new ReadNewickTree(treefile); + int readOk = read->read(ct); - T[0]->assembleTree(); - - /***************************************************/ - // create ouptut tree - respecting pickedGroups // - /***************************************************/ - Tree* outputTree = new Tree(globaldata->Groups.size()); - - outputTree->getSubTree(T[0], globaldata->Groups); - outputTree->assembleTree(); + if (readOk != 0) { m->mothurOut("Read Terminated."); m->mothurOutEndLine(); delete ct; delete read; return 0; } - //no longer need original tree, we have output tree to use and label - for (int i = 0; i < T.size(); i++) { delete T[i]; } globaldata->gTree.clear(); - + vector T = read->getTrees(); + + delete read; + + if (m->control_pressed) { + if (designfile != "") { delete designMap; } + if (sharedfile != "") { for (int i = 0; i < lookup.size(); i++) { delete lookup[i]; } } + else { for (int i = 0; i < lookupFloat.size(); i++) { delete lookupFloat[i]; } } + for (int i = 0; i < T.size(); i++) { delete T[i]; } delete ct; return 0; + } + + T[0]->assembleTree(); + + /***************************************************/ + // create ouptut tree - respecting pickedGroups // + /***************************************************/ + Tree* outputTree = new Tree(m->getNumGroups(), ct); + + outputTree->getSubTree(T[0], m->getGroups()); + outputTree->assembleTree(); - if (m->control_pressed) { - if (sharedfile != "") { for (int i = 0; i < lookup.size(); i++) { delete lookup[i]; } } - else { for (int i = 0; i < lookupFloat.size(); i++) { delete lookupFloat[i]; } } - delete outputTree; delete globaldata->gTreemap; return 0; + //no longer need original tree, we have output tree to use and label + for (int i = 0; i < T.size(); i++) { delete T[i]; } + + if (m->control_pressed) { + if (designfile != "") { delete designMap; } + if (sharedfile != "") { for (int i = 0; i < lookup.size(); i++) { delete lookup[i]; } } + else { for (int i = 0; i < lookupFloat.size(); i++) { delete lookupFloat[i]; } } + delete outputTree; delete ct; return 0; + } + + /***************************************************/ + // get indicator species values // + /***************************************************/ + GetIndicatorSpecies(outputTree); + delete outputTree; delete ct; + + }else { //run with design file only + //get indicator species + GetIndicatorSpecies(); } - /***************************************************/ - // get indicator species values // - /***************************************************/ - GetIndicatorSpecies(outputTree); + if (designfile != "") { delete designMap; } + if (sharedfile != "") { for (int i = 0; i < lookup.size(); i++) { delete lookup[i]; } } + else { for (int i = 0; i < lookupFloat.size(); i++) { delete lookupFloat[i]; } } + + if (m->control_pressed) { for (int i = 0; i < outputNames.size(); i++) { m->mothurRemove(outputNames[i]); } return 0; } - if (m->control_pressed) { - if (sharedfile != "") { for (int i = 0; i < lookup.size(); i++) { delete lookup[i]; } } - else { for (int i = 0; i < lookupFloat.size(); i++) { delete lookupFloat[i]; } } - for (int i = 0; i < outputNames.size(); i++) { remove(outputNames[i].c_str()); } - delete outputTree; delete globaldata->gTreemap; return 0; + //set tree file as new current treefile + if (treefile != "") { + string current = ""; + itTypes = outputTypes.find("tree"); + if (itTypes != outputTypes.end()) { + if ((itTypes->second).size() != 0) { current = (itTypes->second)[0]; m->setTreeFile(current); } + } } + m->mothurOutEndLine(); m->mothurOutEndLine(); m->mothurOut("It took " + toString(time(NULL) - start) + " secs to find the indicator species."); m->mothurOutEndLine(); m->mothurOutEndLine(); m->mothurOut("Output File Names: "); m->mothurOutEndLine(); for (int i = 0; i < outputNames.size(); i++) { m->mothurOut(outputNames[i]); m->mothurOutEndLine(); } @@ -287,6 +399,122 @@ int IndicatorCommand::execute(){ } } //********************************************************************************************************************** +//divide shared or relabund file by groupings in the design file +//report all otu values to file +int IndicatorCommand::GetIndicatorSpecies(){ + try { + string thisOutputDir = outputDir; + if (outputDir == "") { thisOutputDir += m->hasPath(inputFileName); } + map variables; + variables["[filename]"] = thisOutputDir + m->getRootName(m->getSimpleName(inputFileName)); + string outputFileName = getOutputFileName("summary", variables); + outputNames.push_back(outputFileName); outputTypes["summary"].push_back(outputFileName); + + ofstream out; + m->openOutputFile(outputFileName, out); + out.setf(ios::fixed, ios::floatfield); out.setf(ios::showpoint); + m->mothurOutEndLine(); m->mothurOut("Species\tIndicator_Groups\tIndicatorValue\tpValue\n"); + + int numBins = 0; + if (sharedfile != "") { numBins = lookup[0]->getNumBins(); } + else { numBins = lookupFloat[0]->getNumBins(); } + + if (m->control_pressed) { out.close(); return 0; } + + /*****************************************************/ + //create vectors containing rabund info // + /*****************************************************/ + + vector indicatorValues; //size of numBins + vector pValues; + vector indicatorGroups; + map< vector, vector > 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. + + if (sharedfile != "") { + vector< vector > groupings; + set groupsAlreadyAdded; + vector subset; + + //for each grouping + for (int i = 0; i < (designMap->getNamesOfGroups()).size(); i++) { + + for (int k = 0; k < lookup.size(); k++) { + //are you from this grouping? + if (designMap->getGroup(lookup[k]->getGroup()) == (designMap->getNamesOfGroups())[i]) { + subset.push_back(lookup[k]); + groupsAlreadyAdded.insert(lookup[k]->getGroup()); + } + } + if (subset.size() != 0) { groupings.push_back(subset); } + subset.clear(); + } + + if (groupsAlreadyAdded.size() != lookup.size()) { m->mothurOut("[ERROR]: could not make proper groupings."); m->mothurOutEndLine(); } + + indicatorValues = getValues(groupings, indicatorGroups, randomGroupingsMap); + + pValues = getPValues(groupings, lookup.size(), indicatorValues); + }else { + vector< vector > groupings; + set groupsAlreadyAdded; + vector subset; + + //for each grouping + for (int i = 0; i < (designMap->getNamesOfGroups()).size(); i++) { + for (int k = 0; k < lookupFloat.size(); k++) { + //are you from this grouping? + if (designMap->getGroup(lookupFloat[k]->getGroup()) == (designMap->getNamesOfGroups())[i]) { + subset.push_back(lookupFloat[k]); + groupsAlreadyAdded.insert(lookupFloat[k]->getGroup()); + } + } + if (subset.size() != 0) { groupings.push_back(subset); } + subset.clear(); + } + + if (groupsAlreadyAdded.size() != lookupFloat.size()) { m->mothurOut("[ERROR]: could not make proper groupings."); m->mothurOutEndLine(); } + + indicatorValues = getValues(groupings, indicatorGroups, randomGroupingsMap); + + pValues = getPValues(groupings, lookupFloat.size(), indicatorValues); + } + + if (m->control_pressed) { out.close(); return 0; } + + + /******************************************************/ + //output indicator values to table form // + /*****************************************************/ + out << "OTU\tIndicator_Groups\tIndicator_Value\tpValue" << endl; + for (int j = 0; j < indicatorValues.size(); j++) { + + if (m->control_pressed) { out.close(); return 0; } + + out << m->currentBinLabels[j] << '\t' << indicatorGroups[j] << '\t' << indicatorValues[j] << '\t'; + + if (pValues[j] > (1/(float)iters)) { out << pValues[j] << endl; } + else { out << "<" << (1/(float)iters) << endl; } + + if (pValues[j] <= 0.05) { + cout << m->currentBinLabels[j] << '\t' << indicatorGroups[j] << '\t' << indicatorValues[j] << '\t'; + string pValueString = "<" + toString((1/(float)iters)); + if (pValues[j] > (1/(float)iters)) { pValueString = toString(pValues[j]); cout << pValues[j];} + else { cout << "<" << (1/(float)iters); } + m->mothurOutJustToLog(m->currentBinLabels[j] + "\t" + indicatorGroups[j] + "\t" + toString(indicatorValues[j]) + "\t" + pValueString); + m->mothurOutEndLine(); + } + } + + out.close(); + + return 0; + } + catch(exception& e) { + m->errorOut(e, "IndicatorCommand", "GetIndicatorSpecies"); + exit(1); + } +} +//********************************************************************************************************************** //traverse tree finding indicator species values for each otu at each node //label node with otu number that has highest indicator value //report all otu values to file @@ -295,16 +523,30 @@ int IndicatorCommand::GetIndicatorSpecies(Tree*& T){ string thisOutputDir = outputDir; if (outputDir == "") { thisOutputDir += m->hasPath(inputFileName); } - string outputFileName = thisOutputDir + m->getRootName(m->getSimpleName(inputFileName)) + "indicator.summary"; + map variables; + variables["[filename]"] = thisOutputDir + m->getRootName(m->getSimpleName(inputFileName)); + string outputFileName = getOutputFileName("summary",variables); outputNames.push_back(outputFileName); outputTypes["summary"].push_back(outputFileName); ofstream out; m->openOutputFile(outputFileName, out); out.setf(ios::fixed, ios::floatfield); out.setf(ios::showpoint); - + + int numBins = 0; + if (sharedfile != "") { numBins = lookup[0]->getNumBins(); } + else { numBins = lookupFloat[0]->getNumBins(); } + + //print headings + out << "TreeNode\t"; + for (int i = 0; i < numBins; i++) { out << m->currentBinLabels[i] << "_IndGroups" << '\t' << m->currentBinLabels[i] << "_IndValue" << '\t' << "pValue" << '\t'; } + out << endl; + + m->mothurOutEndLine(); m->mothurOut("Node\tSpecies\tIndicator_Groups\tIndicatorValue\tpValue\n"); + string treeOutputDir = outputDir; if (outputDir == "") { treeOutputDir += m->hasPath(treefile); } - string outputTreeFileName = treeOutputDir + m->getRootName(m->getSimpleName(treefile)) + "indicator.tre"; + variables["[filename]"] = treeOutputDir + m->getRootName(m->getSimpleName(treefile)); + string outputTreeFileName = getOutputFileName("tree", variables); //create a map from tree node index to names of descendants, save time later to know which sharedRabund you need @@ -322,7 +564,7 @@ int IndicatorCommand::GetIndicatorSpecies(Tree*& T){ //for each node for (int i = T->getNumLeaves(); i < T->getNumNodes(); i++) { - + //cout << endl << i+1 << endl; if (m->control_pressed) { out.close(); return 0; } /*****************************************************/ @@ -330,6 +572,9 @@ int IndicatorCommand::GetIndicatorSpecies(Tree*& T){ /*****************************************************/ vector indicatorValues; //size of numBins + vector pValues; + vector indicatorGroups; + map< vector, vector > 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. if (sharedfile != "") { vector< vector > groupings; @@ -379,12 +624,10 @@ int IndicatorCommand::GetIndicatorSpecies(Tree*& T){ } if (groupsAlreadyAdded.size() != lookup.size()) { m->mothurOut("[ERROR]: could not make proper groupings."); m->mothurOutEndLine(); } - //for (int k = 0; k < lookup.size(); k++) { - // if (groupsAlreadyAdded.count(lookup[k]->getGroup()) == 0) { cout << lookup[k]->getGroup() << endl; } - //} - - indicatorValues = getValues(groupings); + + indicatorValues = getValues(groupings, indicatorGroups, randomGroupingsMap); + pValues = getPValues(groupings, lookup.size(), indicatorValues); }else { vector< vector > groupings; @@ -431,7 +674,9 @@ int IndicatorCommand::GetIndicatorSpecies(Tree*& T){ if (groupsAlreadyAdded.size() != lookupFloat.size()) { m->mothurOut("[ERROR]: could not make proper groupings."); m->mothurOutEndLine(); } - indicatorValues = getValues(groupings); + indicatorValues = getValues(groupings, indicatorGroups, randomGroupingsMap); + + pValues = getPValues(groupings, lookupFloat.size(), indicatorValues); } if (m->control_pressed) { out.close(); return 0; } @@ -445,12 +690,24 @@ int IndicatorCommand::GetIndicatorSpecies(Tree*& T){ if (m->control_pressed) { out.close(); return 0; } - out << indicatorValues[j] << '\t'; + if (pValues[j] < (1/(float)iters)) { + out << indicatorGroups[j] << '\t' << indicatorValues[j] << '\t' << '<' << (1/(float)iters) << '\t'; + }else { + out << indicatorGroups[j] << '\t' << indicatorValues[j] << '\t' << pValues[j] << '\t'; + } + + if (pValues[j] <= 0.05) { + cout << i+1 << '\t' << m->currentBinLabels[j] << '\t' << indicatorGroups[j] << '\t' << indicatorValues[j] << '\t'; + string pValueString = "<" + toString((1/(float)iters)); + if (pValues[j] > (1/(float)iters)) { pValueString = toString(pValues[j]); cout << pValues[j];} + else { cout << "<" << (1/(float)iters); } + m->mothurOutJustToLog(toString(i) + "\t" + m->currentBinLabels[j] + "\t" + indicatorGroups[j] + "\t" + toString(indicatorValues[j]) + "\t" + pValueString); + m->mothurOutEndLine(); + } } out << endl; T->tree[i].setLabel((i+1)); - } out.close(); @@ -469,10 +726,25 @@ int IndicatorCommand::GetIndicatorSpecies(Tree*& T){ } } //********************************************************************************************************************** -vector IndicatorCommand::getValues(vector< vector >& groupings){ +vector IndicatorCommand::getValues(vector< vector >& groupings, vector& indicatorGroupings, map< vector, vector > groupingsMap){ try { vector values; - + map< vector, vector >::iterator it; + + indicatorGroupings.clear(); + + //create grouping strings + vector groupingsGroups; + for (int j = 0; j < groupings.size(); j++) { + string tempGrouping = ""; + for (int k = 0; k < groupings[j].size()-1; k++) { + tempGrouping += groupings[j][k]->getGroup() + "-"; + } + tempGrouping += groupings[j][groupings[j].size()-1]->getGroup(); + groupingsGroups.push_back(tempGrouping); + } + + //for each otu for (int i = 0; i < groupings[0][0]->getNumBins(); i++) { @@ -481,16 +753,27 @@ vector IndicatorCommand::getValues(vector< vector terms; float AijDenominator = 0.0; vector Bij; + //get overall abundance of each grouping for (int j = 0; j < groupings.size(); j++) { float totalAbund = 0; int numNotZero = 0; for (int k = 0; k < groupings[j].size(); k++) { - totalAbund += groupings[j][k]->getAbundance(i); - if (groupings[j][k]->getAbundance(i) != 0) { numNotZero++; } + vector temp; temp.push_back(j); temp.push_back(k); + it = groupingsMap.find(temp); + + if (it == groupingsMap.end()) { //this one didnt get moved + totalAbund += groupings[j][k]->getAbundance(i); + if (groupings[j][k]->getAbundance(i) != 0.0) { numNotZero++; } + }else { + totalAbund += groupings[(it->second)[0]][(it->second)[1]]->getAbundance(i); + if (groupings[(it->second)[0]][(it->second)[1]]->getAbundance(i) != 0.0) { numNotZero++; } + } + } + //mean abundance float Aij = (totalAbund / (float) groupings[j].size()); terms.push_back(Aij); @@ -501,15 +784,17 @@ vector IndicatorCommand::getValues(vector< vector maxIndVal) { maxIndVal = thisValue; } + if (thisValue > maxIndVal) { maxIndVal = thisValue; maxGrouping = groupingsGroups[j]; } } values.push_back(maxIndVal); + indicatorGroupings.push_back(maxGrouping); } return values; @@ -521,16 +806,25 @@ vector IndicatorCommand::getValues(vector< vector IndicatorCommand::getValues(vector< vector >& groupings){ +vector IndicatorCommand::getValues(vector< vector >& groupings, vector& indicatorGroupings, map< vector, vector > groupingsMap){ try { vector values; + map< vector, vector >::iterator it; + + indicatorGroupings.clear(); + + //create grouping strings + vector groupingsGroups; + for (int j = 0; j < groupings.size(); j++) { + string tempGrouping = ""; + for (int k = 0; k < groupings[j].size()-1; k++) { + tempGrouping += groupings[j][k]->getGroup() + "-"; + } + tempGrouping += groupings[j][groupings[j].size()-1]->getGroup(); + groupingsGroups.push_back(tempGrouping); + } + - /*for (int j = 0; j < groupings.size(); j++) { - cout << "grouping " << j << endl; - for (int k = 0; k < groupings[j].size(); k++) { - cout << groupings[j][k]->getGroup() << endl; - } - }*/ //for each otu for (int i = 0; i < groupings[0][0]->getNumBins(); i++) { vector terms; @@ -542,11 +836,20 @@ vector IndicatorCommand::getValues(vector< vector >& int totalAbund = 0.0; int numNotZero = 0; for (int k = 0; k < groupings[j].size(); k++) { - totalAbund += groupings[j][k]->getAbundance(i); - if (groupings[j][k]->getAbundance(i) != 0.0) { numNotZero++; } + vector temp; temp.push_back(j); temp.push_back(k); + it = groupingsMap.find(temp); - } + if (it == groupingsMap.end()) { //this one didnt get moved + totalAbund += groupings[j][k]->getAbundance(i); + if (groupings[j][k]->getAbundance(i) != 0.0) { numNotZero++; } + }else { + totalAbund += groupings[(it->second)[0]][(it->second)[1]]->getAbundance(i); + if (groupings[(it->second)[0]][(it->second)[1]]->getAbundance(i) != 0.0) { numNotZero++; } + } + } + + //mean abundance float Aij = (totalAbund / (float) groupings[j].size()); terms.push_back(Aij); @@ -557,15 +860,17 @@ vector IndicatorCommand::getValues(vector< vector >& } float maxIndVal = 0.0; + string maxGrouping = ""; for (int j = 0; j < terms.size(); j++) { - float thisAij = (terms[j] / AijDenominator); + float thisAij = (terms[j] / AijDenominator); //relative abundance float thisValue = thisAij * Bij[j] * 100.0; //save largest - if (thisValue > maxIndVal) { maxIndVal = thisValue; } + if (thisValue > maxIndVal) { maxIndVal = thisValue; maxGrouping = groupingsGroups[j]; } } values.push_back(maxIndVal); + indicatorGroupings.push_back(maxGrouping); } return values; @@ -651,8 +956,18 @@ set IndicatorCommand::getDescendantList(Tree*& T, int i, map temp; temp.insert(i); - names.insert(T->tree[i].getName()); nodes[i] = temp; + + if (designfile == "") { + names.insert(T->tree[i].getName()); + }else { + vector myGroup; myGroup.push_back(T->tree[i].getName()); + vector myReps = designMap->getNamesSeqs(myGroup); + for (int k = 0; k < myReps.size(); k++) { + names.insert(myReps[k]); + } + } + }else{ //your descedants are the combination of your childrens descendants names = descendants[lc]; nodes[i] = nodes[lc]; @@ -825,7 +1140,429 @@ int IndicatorCommand::getSharedFloat(){ m->errorOut(e, "IndicatorCommand", "getShared"); exit(1); } -} +} +//********************************************************************************************************************** +vector IndicatorCommand::driver(vector< vector >& groupings, int num, vector indicatorValues, int numIters){ + try { + vector pvalues; + pvalues.resize(indicatorValues.size(), 0); + vector 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. + + for(int i=0;icontrol_pressed) { break; } + map< vector, vector > groupingsMap = randomizeGroupings(groupings, num); + vector randomIndicatorValues = getValues(groupings, notUsedGroupings, groupingsMap); + + for (int j = 0; j < indicatorValues.size(); j++) { + if (randomIndicatorValues[j] >= indicatorValues[j]) { pvalues[j]++; } + } + } + + return pvalues; + + }catch(exception& e) { + m->errorOut(e, "IndicatorCommand", "driver"); + exit(1); + } +} +//********************************************************************************************************************** +vector IndicatorCommand::getPValues(vector< vector >& groupings, int num, vector indicatorValues){ + try { + vector pvalues; + + if(processors == 1){ + pvalues = driver(groupings, num, indicatorValues, iters); + for (int i = 0; i < pvalues.size(); i++) { pvalues[i] /= (double)iters; } + }else{ + //divide iters between processors + vector procIters; + int numItersPerProcessor = iters / processors; + + //divide iters between processes + for (int h = 0; h < processors; h++) { + if(h == processors - 1){ numItersPerProcessor = iters - h * numItersPerProcessor; } + procIters.push_back(numItersPerProcessor); + } + + vector processIDS; + int process = 1; + +#if defined (__APPLE__) || (__MACH__) || (linux) || (__linux) || (__linux__) || (__unix__) || (__unix) + + //loop through and create all the processes you want + while (process != processors) { + int pid = fork(); + + if (pid > 0) { + processIDS.push_back(pid); //create map from line number to pid so you can append files in correct order later + process++; + }else if (pid == 0){ + pvalues = driver(groupings, num, indicatorValues, procIters[process]); + + //pass pvalues to parent + ofstream out; + string tempFile = toString(getpid()) + ".pvalues.temp"; + m->openOutputFile(tempFile, out); + + //pass values + for (int i = 0; i < pvalues.size(); i++) { + out << pvalues[i] << '\t'; + } + out << endl; + + out.close(); + + exit(0); + }else { + m->mothurOut("[ERROR]: unable to spawn the necessary processes."); m->mothurOutEndLine(); + for (int i = 0; i < processIDS.size(); i++) { kill (processIDS[i], SIGINT); } + exit(0); + } + } + + //do my part + pvalues = driver(groupings, num, indicatorValues, procIters[0]); + + //force parent to wait until all the processes are done + for (int i=0;iopenInputFile(tempFile, in); + + ////// to do /////////// + int numTemp; numTemp = 0; + for (int j = 0; j < pvalues.size(); j++) { + in >> numTemp; m->gobble(in); + pvalues[j] += numTemp; + } + in.close(); m->mothurRemove(tempFile); + } + for (int i = 0; i < pvalues.size(); i++) { pvalues[i] /= (double)iters; } +#else + + //fill in functions + vector pDataArray; + DWORD dwThreadIdArray[processors-1]; + HANDLE hThreadArray[processors-1]; + + //Create processor worker threads. + for( int i=1; i > newGroupings; + + for (int k = 0; k < groupings.size(); k++) { + vector newLookup; + for (int l = 0; l < groupings[k].size(); l++) { + SharedRAbundFloatVector* temp = new SharedRAbundFloatVector(); + temp->setLabel(groupings[k][l]->getLabel()); + temp->setGroup(groupings[k][l]->getGroup()); + newLookup.push_back(temp); + } + newGroupings.push_back(newLookup); + } + + //for each bin + for (int l = 0; l < groupings.size(); l++) { + for (int k = 0; k < groupings[l][0]->getNumBins(); k++) { + 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; } + + for (int j = 0; j < groupings[l].size(); j++) { newGroupings[l][j]->push_back(groupings[l][j]->getAbundance(k), groupings[l][j]->getGroup()); } + } + } + + vector copyIValues = indicatorValues; + + indicatorData* temp = new indicatorData(m, procIters[i], newGroupings, num, copyIValues); + pDataArray.push_back(temp); + processIDS.push_back(i); + + hThreadArray[i-1] = CreateThread(NULL, 0, MyIndicatorThreadFunction, pDataArray[i-1], 0, &dwThreadIdArray[i-1]); + } + + //do my part + pvalues = driver(groupings, num, indicatorValues, procIters[0]); + + //Wait until all threads have terminated. + WaitForMultipleObjects(processors-1, hThreadArray, TRUE, INFINITE); + + //Close all thread handles and free memory allocations. + for(int i=0; i < pDataArray.size(); i++){ + for (int j = 0; j < pDataArray[i]->pvalues.size(); j++) { pvalues[j] += pDataArray[i]->pvalues[j]; } + + for (int l = 0; l < pDataArray[i]->groupings.size(); l++) { + for (int j = 0; j < pDataArray[i]->groupings[l].size(); j++) { delete pDataArray[i]->groupings[l][j]; } + } + + CloseHandle(hThreadArray[i]); + delete pDataArray[i]; + } + + for (int i = 0; i < pvalues.size(); i++) { pvalues[i] /= (double)iters; } +#endif + } + + + return pvalues; + } + catch(exception& e) { + m->errorOut(e, "IndicatorCommand", "getPValues"); + exit(1); + } +} + +//********************************************************************************************************************** +//same as above, just data type difference +vector IndicatorCommand::driver(vector< vector >& groupings, int num, vector indicatorValues, int numIters){ + try { + vector pvalues; + pvalues.resize(indicatorValues.size(), 0); + vector 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. + + for(int i=0;icontrol_pressed) { break; } + map< vector, vector > groupingsMap = randomizeGroupings(groupings, num); + vector randomIndicatorValues = getValues(groupings, notUsedGroupings, groupingsMap); + + for (int j = 0; j < indicatorValues.size(); j++) { + if (randomIndicatorValues[j] >= indicatorValues[j]) { pvalues[j]++; } + } + } + + return pvalues; + + }catch(exception& e) { + m->errorOut(e, "IndicatorCommand", "driver"); + exit(1); + } +} +//********************************************************************************************************************** +//same as above, just data type difference +vector IndicatorCommand::getPValues(vector< vector >& groupings, int num, vector indicatorValues){ + try { + vector pvalues; + + if(processors == 1){ + pvalues = driver(groupings, num, indicatorValues, iters); + for (int i = 0; i < pvalues.size(); i++) { pvalues[i] /= (double)iters; } + }else{ + //divide iters between processors + vector procIters; + int numItersPerProcessor = iters / processors; + + //divide iters between processes + for (int h = 0; h < processors; h++) { + if(h == processors - 1){ numItersPerProcessor = iters - h * numItersPerProcessor; } + procIters.push_back(numItersPerProcessor); + } + + vector processIDS; + int process = 1; + +#if defined (__APPLE__) || (__MACH__) || (linux) || (__linux) || (__linux__) || (__unix__) || (__unix) + + //loop through and create all the processes you want + while (process != processors) { + int pid = fork(); + + if (pid > 0) { + processIDS.push_back(pid); //create map from line number to pid so you can append files in correct order later + process++; + }else if (pid == 0){ + pvalues = driver(groupings, num, indicatorValues, procIters[process]); + + //pass pvalues to parent + ofstream out; + string tempFile = toString(getpid()) + ".pvalues.temp"; + m->openOutputFile(tempFile, out); + + //pass values + for (int i = 0; i < pvalues.size(); i++) { + out << pvalues[i] << '\t'; + } + out << endl; + + out.close(); + + exit(0); + }else { + m->mothurOut("[ERROR]: unable to spawn the necessary processes."); m->mothurOutEndLine(); + for (int i = 0; i < processIDS.size(); i++) { kill (processIDS[i], SIGINT); } + exit(0); + } + } + + //do my part + pvalues = driver(groupings, num, indicatorValues, procIters[0]); + + //force parent to wait until all the processes are done + for (int i=0;iopenInputFile(tempFile, in); + + ////// to do /////////// + int numTemp; numTemp = 0; + for (int j = 0; j < pvalues.size(); j++) { + in >> numTemp; m->gobble(in); + pvalues[j] += numTemp; + } + in.close(); m->mothurRemove(tempFile); + } + for (int i = 0; i < pvalues.size(); i++) { pvalues[i] /= (double)iters; } +#else + + //fill in functions + vector pDataArray; + DWORD dwThreadIdArray[processors-1]; + HANDLE hThreadArray[processors-1]; + + //Create processor worker threads. + for( int i=1; i > newGroupings; + + for (int k = 0; k < groupings.size(); k++) { + vector newLookup; + for (int l = 0; l < groupings[k].size(); l++) { + SharedRAbundFloatVector* temp = new SharedRAbundFloatVector(); + temp->setLabel(groupings[k][l]->getLabel()); + temp->setGroup(groupings[k][l]->getGroup()); + newLookup.push_back(temp); + } + newGroupings.push_back(newLookup); + } + + //for each bin + for (int l = 0; l < groupings.size(); l++) { + for (int k = 0; k < groupings[l][0]->getNumBins(); k++) { + 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; } + + for (int j = 0; j < groupings[l].size(); j++) { newGroupings[l][j]->push_back((float)(groupings[l][j]->getAbundance(k)), groupings[l][j]->getGroup()); } + } + } + + vector copyIValues = indicatorValues; + + indicatorData* temp = new indicatorData(m, procIters[i], newGroupings, num, copyIValues); + pDataArray.push_back(temp); + processIDS.push_back(i); + + hThreadArray[i-1] = CreateThread(NULL, 0, MyIndicatorThreadFunction, pDataArray[i-1], 0, &dwThreadIdArray[i-1]); + } + + //do my part + pvalues = driver(groupings, num, indicatorValues, procIters[0]); + + //Wait until all threads have terminated. + WaitForMultipleObjects(processors-1, hThreadArray, TRUE, INFINITE); + + //Close all thread handles and free memory allocations. + for(int i=0; i < pDataArray.size(); i++){ + for (int j = 0; j < pDataArray[i]->pvalues.size(); j++) { pvalues[j] += pDataArray[i]->pvalues[j]; } + + for (int l = 0; l < pDataArray[i]->groupings.size(); l++) { + for (int j = 0; j < pDataArray[i]->groupings[l].size(); j++) { delete pDataArray[i]->groupings[l][j]; } + } + + CloseHandle(hThreadArray[i]); + delete pDataArray[i]; + } + + for (int i = 0; i < pvalues.size(); i++) { pvalues[i] /= (double)iters; } +#endif + } + + + return pvalues; + } + catch(exception& e) { + m->errorOut(e, "IndicatorCommand", "getPValues"); + exit(1); + } +} +//********************************************************************************************************************** +//swap groups between groupings, in essence randomizing the second column of the design file +map< vector, vector > IndicatorCommand::randomizeGroupings(vector< vector >& groupings, int numLookupGroups){ + try { + + map< vector, vector > randomGroupings; + + for (int i = 0; i < numLookupGroups; i++) { + if (m->control_pressed) {break;} + + //get random groups to swap to switch with + //generate random int between 0 and groupings.size()-1 + int z = m->getRandomIndex(groupings.size()-1); + int x = m->getRandomIndex(groupings.size()-1); + int a = m->getRandomIndex(groupings[z].size()-1); + int b = m->getRandomIndex(groupings[x].size()-1); + //cout << i << '\t' << z << '\t' << x << '\t' << a << '\t' << b << endl; + //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; } + + vector from; + vector to; + + from.push_back(z); from.push_back(a); + to.push_back(x); to.push_back(b); + + randomGroupings[from] = to; + } + //cout << "done" << endl; + return randomGroupings; + } + catch(exception& e) { + m->errorOut(e, "IndicatorCommand", "randomizeGroupings"); + exit(1); + } +} +//********************************************************************************************************************** +//swap groups between groupings, in essence randomizing the second column of the design file +map< vector, vector > IndicatorCommand::randomizeGroupings(vector< vector >& groupings, int numLookupGroups){ + try { + + map< vector, vector > randomGroupings; + + for (int i = 0; i < numLookupGroups; i++) { + + //get random groups to swap to switch with + //generate random int between 0 and groupings.size()-1 + int z = m->getRandomIndex(groupings.size()-1); + int x = m->getRandomIndex(groupings.size()-1); + int a = m->getRandomIndex(groupings[z].size()-1); + int b = m->getRandomIndex(groupings[x].size()-1); + //cout << i << '\t' << z << '\t' << x << '\t' << a << '\t' << b << endl; + + vector from; + vector to; + + from.push_back(z); from.push_back(a); + to.push_back(x); to.push_back(b); + + randomGroupings[from] = to; + } + + return randomGroupings; + } + catch(exception& e) { + m->errorOut(e, "IndicatorCommand", "randomizeGroupings"); + exit(1); + } +} + /*****************************************************************/