unifracweightedcommand.cpp

   1 /*
   2  *  unifracweightedcommand.cpp
   3  *  Mothur
   4  *
   5  *  Created by Sarah Westcott on 2/9/09.
   6  *  Copyright 2009 Schloss Lab UMASS Amherst. All rights reserved.
   7  *
   8  */
   9
  10 #include "unifracweightedcommand.h"
  11 #include "consensus.h"
  12 #include "subsample.h"
  13 #include "treereader.h"
  14
  15 //**********************************************************************************************************************
  16 vector<string> UnifracWeightedCommand::setParameters(){
  17         try {
  18                 CommandParameter ptree("tree", "InputTypes", "", "", "none", "none", "none","weighted-wsummary",false,true,true); parameters.push_back(ptree);
  19         CommandParameter pname("name", "InputTypes", "", "", "NameCount", "none", "none","",false,false,true); parameters.push_back(pname);
  20         CommandParameter pcount("count", "InputTypes", "", "", "NameCount-CountGroup", "none", "none","",false,false,true); parameters.push_back(pcount);
  21                 CommandParameter pgroup("group", "InputTypes", "", "", "CountGroup", "none", "none","",false,false,true); parameters.push_back(pgroup);
  22                 CommandParameter pgroups("groups", "String", "", "", "", "", "","",false,false); parameters.push_back(pgroups);
  23                 CommandParameter piters("iters", "Number", "", "1000", "", "", "","",false,false); parameters.push_back(piters);
  24                 CommandParameter pprocessors("processors", "Number", "", "1", "", "", "","",false,false,true); parameters.push_back(pprocessors);
  25         CommandParameter psubsample("subsample", "String", "", "", "", "", "","",false,false); parameters.push_back(psubsample);
  26         CommandParameter pconsensus("consensus", "Boolean", "", "F", "", "", "","tree",false,false); parameters.push_back(pconsensus);
  27         CommandParameter prandom("random", "Boolean", "", "F", "", "", "","",false,false); parameters.push_back(prandom);
  28                 CommandParameter pdistance("distance", "Multiple", "column-lt-square-phylip", "column", "", "", "","phylip-column",false,false); parameters.push_back(pdistance);
  29                 CommandParameter proot("root", "Boolean", "F", "", "", "", "","",false,false); parameters.push_back(proot);
  30                 CommandParameter pinputdir("inputdir", "String", "", "", "", "", "","",false,false); parameters.push_back(pinputdir);
  31                 CommandParameter poutputdir("outputdir", "String", "", "", "", "", "","",false,false); parameters.push_back(poutputdir);
  32
  33                 vector<string> myArray;
  34                 for (int i = 0; i < parameters.size(); i++) {   myArray.push_back(parameters[i].name);          }
  35                 return myArray;
  36         }
  37         catch(exception& e) {
  38                 m->errorOut(e, "UnifracWeightedCommand", "setParameters");
  39                 exit(1);
  40         }
  41 }
  42 //**********************************************************************************************************************
  43 string UnifracWeightedCommand::getHelpString(){
  44         try {
  45                 string helpString = "";
  46                 helpString += "The unifrac.weighted command parameters are tree, group, name, count, groups, iters, distance, processors, root, subsample, consensus and random.  tree parameter is required unless you have valid current tree file.\n";
  47                 helpString += "The groups parameter allows you to specify which of the groups in your groupfile you would like analyzed.  You must enter at least 2 valid groups.\n";
  48                 helpString += "The group names are separated by dashes.  The iters parameter allows you to specify how many random trees you would like compared to your tree.\n";
  49                 helpString += "The distance parameter allows you to create a distance file from the results. The default is false.\n";
  50                 helpString += "The random parameter allows you to shut off the comparison to random trees. The default is false, meaning don't compare your trees with randomly generated trees.\n";
  51                 helpString += "The root parameter allows you to include the entire root in your calculations. The default is false, meaning stop at the root for this comparision instead of the root of the entire tree.\n";
  52                 helpString += "The processors parameter allows you to specify the number of processors to use. The default is 1.\n";
  53         helpString += "The subsample parameter allows you to enter the size pergroup of the sample or you can set subsample=T and mothur will use the size of your smallest group. The subsample parameter may only be used with a group file.\n";
  54         helpString += "The consensus parameter allows you to indicate you would like trees built from distance matrices created with the results, as well as a consensus tree built from these trees. Default=F.\n";
  55         helpString += "The unifrac.weighted command should be in the following format: unifrac.weighted(groups=yourGroups, iters=yourIters).\n";
  56                 helpString += "Example unifrac.weighted(groups=A-B-C, iters=500).\n";
  57                 helpString += "The default value for groups is all the groups in your groupfile, and iters is 1000.\n";
  58                 helpString += "The unifrac.weighted command output two files: .weighted and .wsummary their descriptions are in the manual.\n";
  59                 helpString += "Note: No spaces between parameter labels (i.e. groups), '=' and parameters (i.e.yourGroups).\n";
  60                 return helpString;
  61         }
  62         catch(exception& e) {
  63                 m->errorOut(e, "UnifracWeightedCommand", "getHelpString");
  64                 exit(1);
  65         }
  66 }
  67 //**********************************************************************************************************************
  68 string UnifracWeightedCommand::getOutputPattern(string type) {
  69     try {
  70         string pattern = "";
  71         if (type == "weighted")            {  pattern = "[filename],weighted-[filename],[tag],weighted";   }
  72         else if (type == "wsummary")        {  pattern = "[filename],wsummary";   }
  73         else if (type == "phylip")           {  pattern = "[filename],[tag],[tag2],dist";   }
  74         else if (type == "column")           {  pattern = "[filename],[tag],[tag2],dist";   }
  75         else if (type == "tree")             {  pattern = "[filename],[tag],[tag2],tre";   }
  76         else { m->mothurOut("[ERROR]: No definition for type " + type + " output pattern.\n"); m->control_pressed = true;  }
  77
  78         return pattern;
  79     }
  80     catch(exception& e) {
  81         m->errorOut(e, "UnifracWeightedCommand", "getOutputPattern");
  82         exit(1);
  83     }
  84 }
  85 //**********************************************************************************************************************
  86 UnifracWeightedCommand::UnifracWeightedCommand(){
  87         try {
  88                 abort = true; calledHelp = true;
  89                 setParameters();
  90                 vector<string> tempOutNames;
  91                 outputTypes["weighted"] = tempOutNames;
  92                 outputTypes["wsummary"] = tempOutNames;
  93                 outputTypes["phylip"] = tempOutNames;
  94                 outputTypes["column"] = tempOutNames;
  95         outputTypes["tree"] = tempOutNames;
  96         }
  97         catch(exception& e) {
  98                 m->errorOut(e, "UnifracWeightedCommand", "UnifracWeightedCommand");
  99                 exit(1);
 100         }
 101 }
 102
 103 /***********************************************************/
 104 UnifracWeightedCommand::UnifracWeightedCommand(string option) {
 105         try {
 106                 abort = false; calledHelp = false;
 107
 108                 //allow user to run help
 109                 if(option == "help") { help(); abort = true; calledHelp = true; }
 110                 else if(option == "citation") { citation(); abort = true; calledHelp = true;}
 111
 112                 else {
 113                         vector<string> myArray = setParameters();
 114
 115                         OptionParser parser(option);
 116                         map<string,string> parameters=parser.getParameters();
 117                         map<string,string>::iterator it;
 118
 119                         ValidParameters validParameter;
 120
 121                         //check to make sure all parameters are valid for command
 122                         for (map<string,string>::iterator it = parameters.begin(); it != parameters.end(); it++) {
 123                                 if (validParameter.isValidParameter(it->first, myArray, it->second) != true) {  abort = true;  }
 124                         }
 125
 126                         //initialize outputTypes
 127                         vector<string> tempOutNames;
 128                         outputTypes["weighted"] = tempOutNames;
 129                         outputTypes["wsummary"] = tempOutNames;
 130                         outputTypes["phylip"] = tempOutNames;
 131                         outputTypes["column"] = tempOutNames;
 132             outputTypes["tree"] = tempOutNames;
 133
 134                         //if the user changes the input directory command factory will send this info to us in the output parameter
 135                         string inputDir = validParameter.validFile(parameters, "inputdir", false);
 136                         if (inputDir == "not found"){   inputDir = "";          }
 137                         else {
 138                                 string path;
 139                                 it = parameters.find("tree");
 140                                 //user has given a template file
 141                                 if(it != parameters.end()){
 142                                         path = m->hasPath(it->second);
 143                                         //if the user has not given a path then, add inputdir. else leave path alone.
 144                                         if (path == "") {       parameters["tree"] = inputDir + it->second;             }
 145                                 }
 146
 147                                 it = parameters.find("group");
 148                                 //user has given a template file
 149                                 if(it != parameters.end()){
 150                                         path = m->hasPath(it->second);
 151                                         //if the user has not given a path then, add inputdir. else leave path alone.
 152                                         if (path == "") {       parameters["group"] = inputDir + it->second;            }
 153                                 }
 154
 155                                 it = parameters.find("name");
 156                                 //user has given a template file
 157                                 if(it != parameters.end()){
 158                                         path = m->hasPath(it->second);
 159                                         //if the user has not given a path then, add inputdir. else leave path alone.
 160                                         if (path == "") {       parameters["name"] = inputDir + it->second;             }
 161                                 }
 162
 163                 it = parameters.find("count");
 164                                 //user has given a template file
 165                                 if(it != parameters.end()){
 166                                         path = m->hasPath(it->second);
 167                                         //if the user has not given a path then, add inputdir. else leave path alone.
 168                                         if (path == "") {       parameters["count"] = inputDir + it->second;            }
 169                                 }
 170                         }
 171
 172                         //check for required parameters
 173                         treefile = validParameter.validFile(parameters, "tree", true);
 174                         if (treefile == "not open") { treefile = ""; abort = true; }
 175                         else if (treefile == "not found") {                             //if there is a current design file, use it
 176                                 treefile = m->getTreeFile();
 177                                 if (treefile != "") { m->mothurOut("Using " + treefile + " as input file for the tree parameter."); m->mothurOutEndLine(); }
 178                                 else {  m->mothurOut("You have no current tree file and the tree parameter is required."); m->mothurOutEndLine(); abort = true; }
 179                         }else { m->setTreeFile(treefile); }
 180
 181                         //check for required parameters
 182                         groupfile = validParameter.validFile(parameters, "group", true);
 183                         if (groupfile == "not open") { abort = true; }
 184                         else if (groupfile == "not found") { groupfile = ""; }
 185                         else { m->setGroupFile(groupfile); }
 186
 187                         namefile = validParameter.validFile(parameters, "name", true);
 188                         if (namefile == "not open") { namefile = ""; abort = true; }
 189                         else if (namefile == "not found") { namefile = ""; }
 190                         else { m->setNameFile(namefile); }
 191
 192             countfile = validParameter.validFile(parameters, "count", true);
 193                         if (countfile == "not open") { countfile = ""; abort = true; }
 194                         else if (countfile == "not found") { countfile = "";  }
 195                         else { m->setCountTableFile(countfile); }
 196
 197             if ((namefile != "") && (countfile != "")) {
 198                 m->mothurOut("[ERROR]: you may only use one of the following: name or count."); m->mothurOutEndLine(); abort = true;
 199             }
 200
 201             if ((groupfile != "") && (countfile != "")) {
 202                 m->mothurOut("[ERROR]: you may only use one of the following: group or count."); m->mothurOutEndLine(); abort=true;
 203             }
 204
 205                         outputDir = validParameter.validFile(parameters, "outputdir", false);           if (outputDir == "not found"){  outputDir = m->hasPath(treefile);       }
 206
 207
 208                         //check for optional parameter and set defaults
 209                         // ...at some point should added some additional type checking...
 210                         groups = validParameter.validFile(parameters, "groups", false);
 211                         if (groups == "not found") { groups = ""; }
 212                         else {
 213                                 m->splitAtDash(groups, Groups);
 214                                 m->setGroups(Groups);
 215                         }
 216
 217                         itersString = validParameter.validFile(parameters, "iters", false);                     if (itersString == "not found") { itersString = "1000"; }
 218                         m->mothurConvert(itersString, iters);
 219
 220                         string temp = validParameter.validFile(parameters, "distance", false);
 221                         if (temp == "not found") { phylip = false; outputForm = ""; }
 222                         else{
 223                 if (temp=="phylip") { temp = "lt"; }
 224                                 if ((temp == "lt") || (temp == "column") || (temp == "square")) {  phylip = true;  outputForm = temp; }
 225                                 else { m->mothurOut("Options for distance are: lt, square, or column. Using lt."); m->mothurOutEndLine(); phylip = true; outputForm = "lt"; }
 226                         }
 227
 228                         temp = validParameter.validFile(parameters, "random", false);                           if (temp == "not found") { temp = "F"; }
 229                         random = m->isTrue(temp);
 230
 231                         temp = validParameter.validFile(parameters, "root", false);                                     if (temp == "not found") { temp = "F"; }
 232                         includeRoot = m->isTrue(temp);
 233
 234                         temp = validParameter.validFile(parameters, "processors", false);       if (temp == "not found"){       temp = m->getProcessors();      }
 235                         m->setProcessors(temp);
 236                         m->mothurConvert(temp, processors);
 237
 238             temp = validParameter.validFile(parameters, "subsample", false);            if (temp == "not found") { temp = "F"; }
 239                         if (m->isNumeric1(temp)) { m->mothurConvert(temp, subsampleSize); subsample = true; }
 240             else {
 241                 if (m->isTrue(temp)) { subsample = true; subsampleSize = -1; }  //we will set it to smallest group later
 242                 else { subsample = false; }
 243             }
 244
 245             if (!subsample) { subsampleIters = 0;   }
 246             else { subsampleIters = iters;          }
 247
 248             temp = validParameter.validFile(parameters, "consensus", false);                                    if (temp == "not found") { temp = "F"; }
 249                         consensus = m->isTrue(temp);
 250
 251                         if (subsample && random) {  m->mothurOut("[ERROR]: random must be false, if subsample=t.\n"); abort=true;  }
 252                         if (countfile == "") { if (subsample && (groupfile == "")) {  m->mothurOut("[ERROR]: if subsample=t, a group file must be provided.\n"); abort=true;  } }
 253             else {
 254                 CountTable testCt;
 255                 if ((!testCt.testGroups(countfile)) && (subsample)) {
 256                     m->mothurOut("[ERROR]: if subsample=t, a count file with group info must be provided.\n"); abort=true;
 257                 }
 258             }
 259             if (subsample && (!phylip)) { phylip=true; outputForm = "lt"; }
 260             if (consensus && (!subsample)) { m->mothurOut("[ERROR]: you cannot use consensus without subsample.\n"); abort=true; }
 261
 262                         if (countfile=="") {
 263                 if (namefile == "") {
 264                     vector<string> files; files.push_back(treefile);
 265                     parser.getNameFile(files);
 266                 }
 267             }
 268                 }
 269
 270
 271         }
 272         catch(exception& e) {
 273                 m->errorOut(e, "UnifracWeightedCommand", "UnifracWeightedCommand");
 274                 exit(1);
 275         }
 276 }
 277 /***********************************************************/
 278 int UnifracWeightedCommand::execute() {
 279         try {
 280
 281                 if (abort == true) { if (calledHelp) { return 0; }  return 2;   }
 282
 283                 m->setTreeFile(treefile);
 284
 285         TreeReader* reader;
 286         if (countfile == "") { reader = new TreeReader(treefile, groupfile, namefile); }
 287         else { reader = new TreeReader(treefile, countfile); }
 288         T = reader->getTrees();
 289         ct = T[0]->getCountTable();
 290         delete reader;
 291
 292         if (m->control_pressed) {  delete ct; for (int i = 0; i < T.size(); i++) { delete T[i]; } return 0; }
 293
 294         map<string, string> variables;
 295                 variables["[filename]"] = outputDir + m->getSimpleName(treefile);
 296                 sumFile = getOutputFileName("wsummary",variables);
 297                 m->openOutputFile(sumFile, outSum);
 298                 outputNames.push_back(sumFile);  outputTypes["wsummary"].push_back(sumFile);
 299
 300         SharedUtil util;
 301                 string s; //to make work with setgroups
 302                 Groups = m->getGroups();
 303                 vector<string> nameGroups = ct->getNamesOfGroups();
 304         if (nameGroups.size() < 2) { m->mothurOut("[ERROR]: You cannot run unifrac.weighted with less than 2 groups, aborting.\n"); delete ct; for (int i = 0; i < T.size(); i++) { delete T[i]; } return 0; }
 305                 util.setGroups(Groups, nameGroups, s, numGroups, "weighted");   //sets the groups the user wants to analyze
 306                 m->setGroups(Groups);
 307
 308         if (m->control_pressed) {  delete ct; for (int i = 0; i < T.size(); i++) { delete T[i]; } return 0; }
 309
 310                 Weighted weighted(includeRoot);
 311
 312                 int start = time(NULL);
 313
 314         //set or check size
 315         if (subsample) {
 316             //user has not set size, set size = smallest samples size
 317             if (subsampleSize == -1) {
 318                 vector<string> temp; temp.push_back(Groups[0]);
 319                 subsampleSize = ct->getGroupCount(Groups[0]); //num in first group
 320                 for (int i = 1; i < Groups.size(); i++) {
 321                     int thisSize = ct->getGroupCount(Groups[i]);
 322                     if (thisSize < subsampleSize) {     subsampleSize = thisSize;       }
 323                 }
 324                 m->mothurOut("\nSetting subsample size to " + toString(subsampleSize) + ".\n\n");
 325             }else { //eliminate any too small groups
 326                 vector<string> newGroups = Groups;
 327                 Groups.clear();
 328                 for (int i = 0; i < newGroups.size(); i++) {
 329                     int thisSize = ct->getGroupCount(newGroups[i]);
 330
 331                     if (thisSize >= subsampleSize) {    Groups.push_back(newGroups[i]); }
 332                     else {   m->mothurOut("You have selected a size that is larger than "+newGroups[i]+" number of sequences, removing "+newGroups[i]+".\n"); }
 333                 }
 334                 m->setGroups(Groups);
 335             }
 336         }
 337
 338         //here in case some groups are removed by subsample
 339         util.getCombos(groupComb, Groups, numComp);
 340
 341         if (numComp < processors) { processors = numComp; }
 342
 343         if (consensus && (numComp < 2)) { m->mothurOut("consensus can only be used with numComparisions greater than 1, setting consensus=f.\n"); consensus=false; }
 344
 345         //get weighted scores for users trees
 346         for (int i = 0; i < T.size(); i++) {
 347
 348             if (m->control_pressed) { delete ct; for (int i = 0; i < T.size(); i++) { delete T[i]; } outSum.close(); for (int i = 0; i < outputNames.size(); i++) {     m->mothurRemove(outputNames[i]);  } return 0; }
 349
 350             counter = 0;
 351             rScores.resize(numComp);  //data[0] = weightedscore AB, data[1] = weightedscore AC...
 352             uScores.resize(numComp);  //data[0] = weightedscore AB, data[1] = weightedscore AC...
 353
 354             vector<double> userData; userData.resize(numComp,0);  //weighted score info for user tree. data[0] = weightedscore AB, data[1] = weightedscore AC...
 355             vector<double> randomData; randomData.resize(numComp,0); //weighted score info for random trees. data[0] = weightedscore AB, data[1] = weightedscore AC...
 356
 357             if (random) {
 358                 variables["[filename]"] = outputDir + m->getSimpleName(treefile);
 359                 variables["[tag]"] = toString(i+1);
 360                 string wFileName = getOutputFileName("weighted", variables);
 361                 output = new ColumnFile(wFileName, itersString);
 362                                 outputNames.push_back(wFileName); outputTypes["wweighted"].push_back(wFileName);
 363             }
 364
 365             userData = weighted.getValues(T[i], processors, outputDir); //userData[0] = weightedscore
 366             if (m->control_pressed) { delete ct; for (int i = 0; i < T.size(); i++) { delete T[i]; } if (random) { delete output; } outSum.close(); for (int i = 0; i < outputNames.size(); i++) {      m->mothurRemove(outputNames[i]);  } return 0; }
 367
 368             //save users score
 369             for (int s=0; s<numComp; s++) {
 370                 //add users score to vector of user scores
 371                 uScores[s].push_back(userData[s]);
 372                 //save users tree score for summary file
 373                 utreeScores.push_back(userData[s]);
 374             }
 375
 376             if (random) {  runRandomCalcs(T[i], userData); }
 377
 378             //clear data
 379             rScores.clear();
 380             uScores.clear();
 381             validScores.clear();
 382
 383             //subsample loop
 384             vector< vector<double> > calcDistsTotals;  //each iter, each groupCombos dists. this will be used to make .dist files
 385             for (int thisIter = 0; thisIter < subsampleIters; thisIter++) { //subsampleIters=0, if subsample=f.
 386                 if (m->control_pressed) { break; }
 387
 388                 //copy to preserve old one - would do this in subsample but memory cleanup becomes messy.
 389                 CountTable* newCt = new CountTable();
 390
 391                 int sampleTime = 0;
 392                 if (m->debug) { sampleTime = time(NULL); }
 393
 394                 //uses method of setting groups to doNotIncludeMe
 395                 SubSample sample;
 396                 Tree* subSampleTree = sample.getSample(T[i], ct, newCt, subsampleSize);
 397
 398                 if (m->debug) { m->mothurOut("[DEBUG]: iter " + toString(thisIter) + " took " + toString(time(NULL) - sampleTime) + " seconds to sample tree.\n"); }
 399
 400                 //call new weighted function
 401                 vector<double> iterData; iterData.resize(numComp,0);
 402                 Weighted thisWeighted(includeRoot);
 403                 iterData = thisWeighted.getValues(subSampleTree, processors, outputDir); //userData[0] = weightedscore
 404
 405                 //save data to make ave dist, std dist
 406                 calcDistsTotals.push_back(iterData);
 407
 408                 delete newCt;
 409                 delete subSampleTree;
 410
 411                 if((thisIter+1) % 100 == 0){    m->mothurOutJustToScreen(toString(thisIter+1)+"\n");    }
 412             }
 413
 414             if (m->control_pressed) { delete ct; for (int i = 0; i < T.size(); i++) { delete T[i]; } if (random) { delete output; } outSum.close(); for (int i = 0; i < outputNames.size(); i++) {      m->mothurRemove(outputNames[i]);  } return 0; }
 415
 416             if (subsample) {  getAverageSTDMatrices(calcDistsTotals, i); }
 417             if (consensus) {  getConsensusTrees(calcDistsTotals, i);  }
 418         }
 419
 420
 421                 if (m->control_pressed) { delete ct; for (int i = 0; i < T.size(); i++) { delete T[i]; } outSum.close(); for (int i = 0; i < outputNames.size(); i++) { m->mothurRemove(outputNames[i]);  } return 0;  }
 422
 423         if (phylip) {   createPhylipFile();             }
 424
 425                 printWSummaryFile();
 426
 427                 //clear out users groups
 428                 m->clearGroups();
 429                 delete ct;
 430                 for (int i = 0; i < T.size(); i++) { delete T[i]; }
 431
 432                 if (m->control_pressed) { for (int i = 0; i < outputNames.size(); i++) {        m->mothurRemove(outputNames[i]);  } return 0; }
 433
 434                 m->mothurOut("It took " + toString(time(NULL) - start) + " secs to run unifrac.weighted."); m->mothurOutEndLine();
 435
 436                 //set phylip file as new current phylipfile
 437                 string current = "";
 438                 itTypes = outputTypes.find("phylip");
 439                 if (itTypes != outputTypes.end()) {
 440                         if ((itTypes->second).size() != 0) { current = (itTypes->second)[0]; m->setPhylipFile(current); }
 441                 }
 442
 443                 //set column file as new current columnfile
 444                 itTypes = outputTypes.find("column");
 445                 if (itTypes != outputTypes.end()) {
 446                         if ((itTypes->second).size() != 0) { current = (itTypes->second)[0]; m->setColumnFile(current); }
 447                 }
 448
 449                 m->mothurOutEndLine();
 450                 m->mothurOut("Output File Names: "); m->mothurOutEndLine();
 451                 for (int i = 0; i < outputNames.size(); i++) {  m->mothurOut(outputNames[i]); m->mothurOutEndLine();    }
 452                 m->mothurOutEndLine();
 453
 454                 return 0;
 455
 456         }
 457         catch(exception& e) {
 458                 m->errorOut(e, "UnifracWeightedCommand", "execute");
 459                 exit(1);
 460         }
 461 }
 462 /**************************************************************************************************/
 463 int UnifracWeightedCommand::getAverageSTDMatrices(vector< vector<double> >& dists, int treeNum) {
 464         try {
 465         //we need to find the average distance and standard deviation for each groups distance
 466
 467         //finds sum
 468         vector<double> averages = m->getAverages(dists);
 469
 470         //find standard deviation
 471         vector<double> stdDev = m->getStandardDeviation(dists, averages);
 472
 473         //make matrix with scores in it
 474         vector< vector<double> > avedists;      //avedists.resize(m->getNumGroups());
 475         for (int i = 0; i < m->getNumGroups(); i++) {
 476             vector<double> temp;
 477             for (int j = 0; j < m->getNumGroups(); j++) { temp.push_back(0.0); }
 478             avedists.push_back(temp);
 479         }
 480
 481         //make matrix with scores in it
 482         vector< vector<double> > stddists;      //stddists.resize(m->getNumGroups());
 483         for (int i = 0; i < m->getNumGroups(); i++) {
 484             vector<double> temp;
 485             for (int j = 0; j < m->getNumGroups(); j++) { temp.push_back(0.0); }
 486             //stddists[i].resize(m->getNumGroups(), 0.0);
 487             stddists.push_back(temp);
 488         }
 489
 490
 491         //flip it so you can print it
 492         int count = 0;
 493         for (int r=0; r<m->getNumGroups(); r++) {
 494             for (int l = 0; l < r; l++) {
 495                 avedists[r][l] = averages[count];
 496                 avedists[l][r] = averages[count];
 497                 stddists[r][l] = stdDev[count];
 498                 stddists[l][r] = stdDev[count];
 499                 count++;
 500             }
 501         }
 502
 503         map<string, string> variables;
 504                 variables["[filename]"] = outputDir + m->getSimpleName(treefile);
 505         variables["[tag]"] = toString(treeNum+1);
 506         variables["[tag2]"] = "weighted.ave";
 507         string aveFileName = getOutputFileName("phylip",variables);
 508         if (outputForm != "column") { outputNames.push_back(aveFileName); outputTypes["phylip"].push_back(aveFileName);  }
 509         else { outputNames.push_back(aveFileName); outputTypes["column"].push_back(aveFileName);  }
 510         ofstream out;
 511         m->openOutputFile(aveFileName, out);
 512
 513         variables["[tag2]"] = "weighted.std";
 514         string stdFileName = getOutputFileName("phylip",variables);
 515         if (outputForm != "column") { outputNames.push_back(stdFileName); outputTypes["phylip"].push_back(stdFileName); }
 516         else { outputNames.push_back(stdFileName); outputTypes["column"].push_back(stdFileName); }
 517         ofstream outStd;
 518         m->openOutputFile(stdFileName, outStd);
 519
 520         if ((outputForm == "lt") || (outputForm == "square")) {
 521             //output numSeqs
 522             out << m->getNumGroups() << endl;
 523             outStd << m->getNumGroups() << endl;
 524         }
 525
 526         //output to file
 527         for (int r=0; r<m->getNumGroups(); r++) {
 528             //output name
 529             string name = (m->getGroups())[r];
 530             if (name.length() < 10) { //pad with spaces to make compatible
 531                 while (name.length() < 10) {  name += " ";  }
 532             }
 533
 534             if (outputForm == "lt") {
 535                 out << name << '\t';
 536                 outStd << name << '\t';
 537
 538                 //output distances
 539                 for (int l = 0; l < r; l++) {   out  << avedists[r][l] << '\t';  outStd  << stddists[r][l] << '\t';}
 540                 out << endl;  outStd << endl;
 541             }else if (outputForm == "square") {
 542                 out << name << '\t';
 543                 outStd << name << '\t';
 544
 545                 //output distances
 546                 for (int l = 0; l < m->getNumGroups(); l++) {   out  << avedists[r][l] << '\t'; outStd  << stddists[r][l] << '\t'; }
 547                 out << endl; outStd << endl;
 548             }else{
 549                 //output distances
 550                 for (int l = 0; l < r; l++) {
 551                     string otherName = (m->getGroups())[l];
 552                     if (otherName.length() < 10) { //pad with spaces to make compatible
 553                         while (otherName.length() < 10) {  otherName += " ";  }
 554                     }
 555
 556                     out  << name << '\t' << otherName << avedists[r][l] << endl;
 557                     outStd  << name << '\t' << otherName << stddists[r][l] << endl;
 558                 }
 559             }
 560         }
 561         out.close();
 562         outStd.close();
 563
 564         return 0;
 565     }
 566         catch(exception& e) {
 567                 m->errorOut(e, "UnifracWeightedCommand", "getAverageSTDMatrices");
 568                 exit(1);
 569         }
 570 }
 571
 572 /**************************************************************************************************/
 573 int UnifracWeightedCommand::getConsensusTrees(vector< vector<double> >& dists, int treeNum) {
 574         try {
 575
 576         //used in tree constructor
 577         m->runParse = false;
 578
 579         ///create treemap class from groupmap for tree class to use
 580         CountTable newCt;
 581         set<string> nameMap;
 582         map<string, string> groupMap;
 583         set<string> gps;
 584         for (int i = 0; i < m->getGroups().size(); i++) {
 585             nameMap.insert(m->getGroups()[i]);
 586             gps.insert(m->getGroups()[i]);
 587             groupMap[m->getGroups()[i]] = m->getGroups()[i];
 588         }
 589         newCt.createTable(nameMap, groupMap, gps);
 590
 591         //clear  old tree names if any
 592         m->Treenames.clear();
 593
 594         //fills globaldatas tree names
 595         m->Treenames = m->getGroups();
 596
 597         vector<Tree*> newTrees = buildTrees(dists, treeNum, newCt); //also creates .all.tre file containing the trees created
 598
 599         if (m->control_pressed) { return 0; }
 600
 601         Consensus con;
 602         Tree* conTree = con.getTree(newTrees);
 603
 604         //create a new filename
 605         map<string, string> variables;
 606                 variables["[filename]"] = outputDir + m->getRootName(m->getSimpleName(treefile));
 607         variables["[tag]"] = toString(treeNum+1);
 608         variables["[tag2]"] = "weighted.cons";
 609         string conFile = getOutputFileName("tree",variables);
 610         outputNames.push_back(conFile); outputTypes["tree"].push_back(conFile);
 611         ofstream outTree;
 612         m->openOutputFile(conFile, outTree);
 613
 614         if (conTree != NULL) { conTree->print(outTree, "boot"); delete conTree; }
 615         outTree.close();
 616
 617         return 0;
 618
 619     }
 620         catch(exception& e) {
 621                 m->errorOut(e, "UnifracWeightedCommand", "getConsensusTrees");
 622                 exit(1);
 623         }
 624 }
 625 /**************************************************************************************************/
 626
 627 vector<Tree*> UnifracWeightedCommand::buildTrees(vector< vector<double> >& dists, int treeNum, CountTable& myct) {
 628         try {
 629
 630         vector<Tree*> trees;
 631
 632         //create a new filename
 633         map<string, string> variables;
 634                 variables["[filename]"] = outputDir + m->getRootName(m->getSimpleName(treefile));
 635         variables["[tag]"] = toString(treeNum+1);
 636         variables["[tag2]"] = "weighted.all";
 637         string outputFile = getOutputFileName("tree",variables);
 638         outputNames.push_back(outputFile); outputTypes["tree"].push_back(outputFile);
 639
 640         ofstream outAll;
 641         m->openOutputFile(outputFile, outAll);
 642
 643
 644         for (int i = 0; i < dists.size(); i++) { //dists[0] are the dists for the first subsampled tree.
 645
 646             if (m->control_pressed) { break; }
 647
 648             //make matrix with scores in it
 649             vector< vector<double> > sims;      sims.resize(m->getNumGroups());
 650             for (int j = 0; j < m->getNumGroups(); j++) {
 651                 sims[j].resize(m->getNumGroups(), 0.0);
 652             }
 653
 654             int count = 0;
 655                         for (int r=0; r<m->getNumGroups(); r++) {
 656                                 for (int l = 0; l < r; l++) {
 657                     double sim = -(dists[i][count]-1.0);
 658                                         sims[r][l] = sim;
 659                                         sims[l][r] = sim;
 660                                         count++;
 661                                 }
 662                         }
 663
 664             //create tree
 665             Tree* tempTree = new Tree(&myct, sims);
 666             tempTree->assembleTree();
 667
 668             trees.push_back(tempTree);
 669
 670             //print tree
 671             tempTree->print(outAll);
 672         }
 673
 674         outAll.close();
 675
 676         if (m->control_pressed) {  for (int i = 0; i < trees.size(); i++) {  delete trees[i]; trees[i] = NULL; } m->mothurRemove(outputFile); }
 677
 678         return trees;
 679     }
 680         catch(exception& e) {
 681                 m->errorOut(e, "UnifracWeightedCommand", "buildTrees");
 682                 exit(1);
 683         }
 684 }
 685 /**************************************************************************************************/
 686
 687 int UnifracWeightedCommand::runRandomCalcs(Tree* thisTree, vector<double> usersScores) {
 688         try {
 689
 690         //calculate number of comparisons i.e. with groups A,B,C = AB, AC, BC = 3;
 691         vector< vector<string> > namesOfGroupCombos;
 692         for (int a=0; a<numGroups; a++) {
 693             for (int l = 0; l < a; l++) {
 694                 vector<string> groups; groups.push_back((m->getGroups())[a]); groups.push_back((m->getGroups())[l]);
 695                 namesOfGroupCombos.push_back(groups);
 696             }
 697         }
 698
 699         lines.clear();
 700
 701         //breakdown work between processors
 702         int numPairs = namesOfGroupCombos.size();
 703         int numPairsPerProcessor = numPairs / processors;
 704
 705         for (int i = 0; i < processors; i++) {
 706             int startPos = i * numPairsPerProcessor;
 707             if(i == processors - 1){ numPairsPerProcessor = numPairs - i * numPairsPerProcessor; }
 708             lines.push_back(linePair(startPos, numPairsPerProcessor));
 709         }
 710
 711
 712         //get scores for random trees
 713         for (int j = 0; j < iters; j++) {
 714 //#if defined (__APPLE__) || (__MACH__) || (linux) || (__linux) || (__linux__) || (__unix__) || (__unix)
 715             //if(processors == 1){
 716               //  driver(thisTree,  namesOfGroupCombos, 0, namesOfGroupCombos.size(),  rScores);
 717            // }else{
 718                 createProcesses(thisTree,  namesOfGroupCombos, rScores);
 719            // }
 720 //#else
 721             //driver(thisTree, namesOfGroupCombos, 0, namesOfGroupCombos.size(), rScores);
 722 //#endif
 723
 724
 725             if (m->control_pressed) { delete ct;  for (int i = 0; i < T.size(); i++) { delete T[i]; } delete output; outSum.close(); for (int i = 0; i < outputNames.size(); i++) {     m->mothurRemove(outputNames[i]);  } return 0; }
 726
 727         }
 728         lines.clear();
 729
 730         //find the signifigance of the score for summary file
 731         for (int f = 0; f < numComp; f++) {
 732             //sort random scores
 733             sort(rScores[f].begin(), rScores[f].end());
 734
 735             //the index of the score higher than yours is returned
 736             //so if you have 1000 random trees the index returned is 100
 737             //then there are 900 trees with a score greater then you.
 738             //giving you a signifigance of 0.900
 739             int index = findIndex(usersScores[f], f);    if (index == -1) { m->mothurOut("error in UnifracWeightedCommand"); m->mothurOutEndLine(); exit(1); } //error code
 740
 741             //the signifigance is the number of trees with the users score or higher
 742             WScoreSig.push_back((iters-index)/(float)iters);
 743         }
 744
 745         //out << "Tree# " << i << endl;
 746         calculateFreqsCumuls();
 747         printWeightedFile();
 748
 749         delete output;
 750
 751         return 0;
 752     }
 753         catch(exception& e) {
 754                 m->errorOut(e, "UnifracWeightedCommand", "runRandomCalcs");
 755                 exit(1);
 756         }
 757 }
 758 /**************************************************************************************************/
 759
 760 int UnifracWeightedCommand::createProcesses(Tree* t, vector< vector<string> > namesOfGroupCombos, vector< vector<double> >& scores) {
 761         try {
 762         int process = 1;
 763                 vector<int> processIDS;
 764                 EstOutput results;
 765
 766 #if defined (__APPLE__) || (__MACH__) || (linux) || (__linux) || (__linux__) || (__unix__) || (__unix)
 767                 //loop through and create all the processes you want
 768                 while (process != processors) {
 769                         int pid = fork();
 770
 771                         if (pid > 0) {
 772                                 processIDS.push_back(pid);  //create map from line number to pid so you can append files in correct order later
 773                                 process++;
 774                         }else if (pid == 0){
 775                                 driver(t, namesOfGroupCombos, lines[process].start, lines[process].num, scores);
 776
 777                                 //pass numSeqs to parent
 778                                 ofstream out;
 779                                 string tempFile = outputDir + toString(getpid()) + ".weightedcommand.results.temp";
 780                                 m->openOutputFile(tempFile, out);
 781                                 for (int i = lines[process].start; i < (lines[process].start + lines[process].num); i++) { out << scores[i][(scores[i].size()-1)] << '\t';  } out << endl;
 782                                 out.close();
 783
 784                                 exit(0);
 785                         }else {
 786                                 m->mothurOut("[ERROR]: unable to spawn the necessary processes."); m->mothurOutEndLine();
 787                                 for (int i = 0; i < processIDS.size(); i++) { kill (processIDS[i], SIGINT); }
 788                                 exit(0);
 789                         }
 790                 }
 791
 792                 driver(t, namesOfGroupCombos, lines[0].start, lines[0].num, scores);
 793
 794                 //force parent to wait until all the processes are done
 795                 for (int i=0;i<(processors-1);i++) {
 796                         int temp = processIDS[i];
 797                         wait(&temp);
 798                 }
 799
 800                 //get data created by processes
 801                 for (int i=0;i<(processors-1);i++) {
 802
 803                         ifstream in;
 804                         string s = outputDir + toString(processIDS[i]) + ".weightedcommand.results.temp";
 805                         m->openInputFile(s, in);
 806
 807                         double tempScore;
 808                         for (int j = lines[(i+1)].start; j < (lines[(i+1)].start + lines[(i+1)].num); j++) { in >> tempScore; scores[j].push_back(tempScore); }
 809                         in.close();
 810                         m->mothurRemove(s);
 811                 }
 812 #else
 813         //fill in functions
 814         vector<weightedRandomData*> pDataArray;
 815                 DWORD   dwThreadIdArray[processors-1];
 816                 HANDLE  hThreadArray[processors-1];
 817         vector<CountTable*> cts;
 818         vector<Tree*> trees;
 819
 820                 //Create processor worker threads.
 821                 for( int i=1; i<processors; i++ ){
 822             CountTable* copyCount = new CountTable();
 823             copyCount->copy(ct);
 824             Tree* copyTree = new Tree(copyCount);
 825             copyTree->getCopy(t);
 826
 827             cts.push_back(copyCount);
 828             trees.push_back(copyTree);
 829
 830             vector< vector<double> > copyScores = rScores;
 831
 832             weightedRandomData* tempweighted = new weightedRandomData(m, lines[i].start, lines[i].num, namesOfGroupCombos, copyTree, copyCount, includeRoot, copyScores);
 833                         pDataArray.push_back(tempweighted);
 834                         processIDS.push_back(i);
 835
 836                         hThreadArray[i-1] = CreateThread(NULL, 0, MyWeightedRandomThreadFunction, pDataArray[i-1], 0, &dwThreadIdArray[i-1]);
 837                 }
 838
 839                 driver(t, namesOfGroupCombos, lines[0].start, lines[0].num, scores);
 840
 841                 //Wait until all threads have terminated.
 842                 WaitForMultipleObjects(processors-1, hThreadArray, TRUE, INFINITE);
 843
 844                 //Close all thread handles and free memory allocations.
 845                 for(int i=0; i < pDataArray.size(); i++){
 846             for (int j = pDataArray[i]->start; j < (pDataArray[i]->start+pDataArray[i]->num); j++) {
 847                 scores[j].push_back(pDataArray[i]->scores[j][pDataArray[i]->scores[j].size()-1]);
 848             }
 849                         delete cts[i];
 850             delete trees[i];
 851                         CloseHandle(hThreadArray[i]);
 852                         delete pDataArray[i];
 853                 }
 854
 855
 856 #endif
 857
 858         return 0;
 859         }
 860         catch(exception& e) {
 861                 m->errorOut(e, "UnifracWeightedCommand", "createProcesses");
 862                 exit(1);
 863         }
 864 }
 865
 866 /**************************************************************************************************/
 867 int UnifracWeightedCommand::driver(Tree* t, vector< vector<string> > namesOfGroupCombos, int start, int num, vector< vector<double> >& scores) {
 868  try {
 869                 Tree* randT = new Tree(ct);
 870
 871         Weighted weighted(includeRoot);
 872
 873                 for (int h = start; h < (start+num); h++) {
 874
 875                         if (m->control_pressed) { return 0; }
 876
 877                         //initialize weighted score
 878                         string groupA = namesOfGroupCombos[h][0];
 879                         string groupB = namesOfGroupCombos[h][1];
 880
 881                         //copy T[i]'s info.
 882                         randT->getCopy(t);
 883
 884                         //create a random tree with same topology as T[i], but different labels
 885                         randT->assembleRandomUnifracTree(groupA, groupB);
 886
 887                         if (m->control_pressed) { delete randT;  return 0;  }
 888
 889                         //get wscore of random tree
 890                         EstOutput randomData = weighted.getValues(randT, groupA, groupB);
 891
 892                         if (m->control_pressed) { delete randT;  return 0;  }
 893
 894                         //save scores
 895                         scores[h].push_back(randomData[0]);
 896                 }
 897
 898                 delete randT;
 899
 900                 return 0;
 901
 902         }
 903         catch(exception& e) {
 904                 m->errorOut(e, "UnifracWeightedCommand", "driver");
 905                 exit(1);
 906         }
 907 }
 908 /***********************************************************/
 909 void UnifracWeightedCommand::printWeightedFile() {
 910         try {
 911                 vector<double> data;
 912                 vector<string> tags;
 913                 tags.push_back("Score"); tags.push_back("RandFreq"); tags.push_back("RandCumul");
 914
 915                 for(int a = 0; a < numComp; a++) {
 916                         output->initFile(groupComb[a], tags);
 917                         //print each line
 918                         for (map<double,double>::iterator it = validScores.begin(); it != validScores.end(); it++) {
 919                                 data.push_back(it->first);  data.push_back(rScoreFreq[a][it->first]); data.push_back(rCumul[a][it->first]);
 920                                 output->output(data);
 921                                 data.clear();
 922                         }
 923                         output->resetFile();
 924                 }
 925         }
 926         catch(exception& e) {
 927                 m->errorOut(e, "UnifracWeightedCommand", "printWeightedFile");
 928                 exit(1);
 929         }
 930 }
 931
 932
 933 /***********************************************************/
 934 void UnifracWeightedCommand::printWSummaryFile() {
 935         try {
 936                 //column headers
 937                 outSum << "Tree#" << '\t' << "Groups" << '\t' << "WScore" << '\t';
 938                 m->mothurOut("Tree#\tGroups\tWScore\t");
 939                 if (random) { outSum << "WSig"; m->mothurOut("WSig"); }
 940                 outSum << endl; m->mothurOutEndLine();
 941
 942                 //format output
 943                 outSum.setf(ios::fixed, ios::floatfield); outSum.setf(ios::showpoint);
 944
 945                 //print each line
 946                 int count = 0;
 947                 for (int i = 0; i < T.size(); i++) {
 948                         for (int j = 0; j < numComp; j++) {
 949                                 if (random) {
 950                                         if (WScoreSig[count] > (1/(float)iters)) {
 951                                                 outSum << setprecision(6) << i+1 << '\t' << groupComb[j] << '\t' << utreeScores[count] << '\t' << setprecision(itersString.length()) << WScoreSig[count] << endl;
 952                                                 cout << setprecision(6) << i+1 << '\t' << groupComb[j] << '\t' << utreeScores[count] << '\t' << setprecision(itersString.length()) << WScoreSig[count] << endl;
 953                                                 m->mothurOutJustToLog(toString(i+1) +"\t" + groupComb[j] +"\t" + toString(utreeScores[count]) +"\t" +  toString(WScoreSig[count]) + "\n");
 954                                         }else{
 955                                                 outSum << setprecision(6) << i+1 << '\t' << groupComb[j] << '\t' << utreeScores[count] << '\t' << setprecision(itersString.length()) << "<" << (1/float(iters)) << endl;
 956                                                 cout << setprecision(6) << i+1 << '\t' << groupComb[j] << '\t' << utreeScores[count] << '\t' << setprecision(itersString.length()) << "<" << (1/float(iters)) << endl;
 957                                                 m->mothurOutJustToLog(toString(i+1) +"\t" + groupComb[j] +"\t" + toString(utreeScores[count]) +"\t<" +  toString((1/float(iters))) + "\n");
 958                                         }
 959                                 }else{
 960                                         outSum << setprecision(6) << i+1 << '\t' << groupComb[j] << '\t' << utreeScores[count] << endl;
 961                                         cout << setprecision(6) << i+1 << '\t' << groupComb[j] << '\t' << utreeScores[count]  << endl;
 962                                         m->mothurOutJustToLog(toString(i+1) +"\t" + groupComb[j] +"\t" + toString(utreeScores[count]) +"\n");
 963                                 }
 964                                 count++;
 965                         }
 966                 }
 967                 outSum.close();
 968         }
 969         catch(exception& e) {
 970                 m->errorOut(e, "UnifracWeightedCommand", "printWSummaryFile");
 971                 exit(1);
 972         }
 973 }
 974 /***********************************************************/
 975 void UnifracWeightedCommand::createPhylipFile() {
 976         try {
 977                 int count = 0;
 978                 //for each tree
 979                 for (int i = 0; i < T.size(); i++) {
 980
 981             string phylipFileName;
 982                         map<string, string> variables;
 983             variables["[filename]"] = outputDir + m->getSimpleName(treefile);
 984             variables["[tag]"] = toString(i+1);
 985             if ((outputForm == "lt") || (outputForm == "square")) {
 986                 variables["[tag2]"] = "weighted.phylip";
 987                 phylipFileName = getOutputFileName("phylip",variables);
 988                 outputNames.push_back(phylipFileName); outputTypes["phylip"].push_back(phylipFileName);
 989             }else { //column
 990                 variables["[tag2]"] = "weighted.column";
 991                 phylipFileName = getOutputFileName("column",variables);
 992                 outputNames.push_back(phylipFileName); outputTypes["column"].push_back(phylipFileName);
 993             }
 994
 995
 996                         ofstream out;
 997                         m->openOutputFile(phylipFileName, out);
 998
 999                         if ((outputForm == "lt") || (outputForm == "square")) {
1000                                 //output numSeqs
1001                                 out << m->getNumGroups() << endl;
1002                         }
1003
1004                         //make matrix with scores in it
1005                         vector< vector<float> > dists;  dists.resize(m->getNumGroups());
1006                         for (int i = 0; i < m->getNumGroups(); i++) {
1007                                 dists[i].resize(m->getNumGroups(), 0.0);
1008                         }
1009
1010                         //flip it so you can print it
1011                         for (int r=0; r<m->getNumGroups(); r++) {
1012                                 for (int l = 0; l < r; l++) {
1013                                         dists[r][l] = utreeScores[count];
1014                                         dists[l][r] = utreeScores[count];
1015                                         count++;
1016                                 }
1017                         }
1018
1019                         //output to file
1020                         for (int r=0; r<m->getNumGroups(); r++) {
1021                                 //output name
1022                                 string name = (m->getGroups())[r];
1023                                 if (name.length() < 10) { //pad with spaces to make compatible
1024                                         while (name.length() < 10) {  name += " ";  }
1025                                 }
1026
1027                                 if (outputForm == "lt") {
1028                                         out << name << '\t';
1029
1030                                         //output distances
1031                                         for (int l = 0; l < r; l++) {   out  << dists[r][l] << '\t';  }
1032                                         out << endl;
1033                                 }else if (outputForm == "square") {
1034                                         out << name << '\t';
1035
1036                                         //output distances
1037                                         for (int l = 0; l < m->getNumGroups(); l++) {   out  << dists[r][l] << '\t';  }
1038                                         out << endl;
1039                                 }else{
1040                                         //output distances
1041                                         for (int l = 0; l < r; l++) {
1042                                                 string otherName = (m->getGroups())[l];
1043                                                 if (otherName.length() < 10) { //pad with spaces to make compatible
1044                                                         while (otherName.length() < 10) {  otherName += " ";  }
1045                                                 }
1046
1047                                                 out  << name << '\t' << otherName << dists[r][l] << endl;
1048                                         }
1049                                 }
1050                         }
1051                         out.close();
1052                 }
1053         }
1054         catch(exception& e) {
1055                 m->errorOut(e, "UnifracWeightedCommand", "createPhylipFile");
1056                 exit(1);
1057         }
1058 }
1059 /***********************************************************/
1060 int UnifracWeightedCommand::findIndex(float score, int index) {
1061         try{
1062                 for (int i = 0; i < rScores[index].size(); i++) {
1063                         if (rScores[index][i] >= score) {       return i;       }
1064                 }
1065                 return rScores[index].size();
1066         }
1067         catch(exception& e) {
1068                 m->errorOut(e, "UnifracWeightedCommand", "findIndex");
1069                 exit(1);
1070         }
1071 }
1072
1073 /***********************************************************/
1074
1075 void UnifracWeightedCommand::calculateFreqsCumuls() {
1076         try {
1077                 //clear out old tree values
1078                 rScoreFreq.clear();
1079                 rScoreFreq.resize(numComp);
1080                 rCumul.clear();
1081                 rCumul.resize(numComp);
1082                 validScores.clear();
1083
1084                 //calculate frequency
1085                 for (int f = 0; f < numComp; f++) {
1086                         for (int i = 0; i < rScores[f].size(); i++) { //looks like 0,0,1,1,1,2,4,7...  you want to make a map that say rScoreFreq[0] = 2, rScoreFreq[1] = 3...
1087                                 validScores[rScores[f][i]] = rScores[f][i];
1088                                 map<double,double>::iterator it = rScoreFreq[f].find(rScores[f][i]);
1089                                 if (it != rScoreFreq[f].end()) {
1090                                         rScoreFreq[f][rScores[f][i]]++;
1091                                 }else{
1092                                         rScoreFreq[f][rScores[f][i]] = 1;
1093                                 }
1094                         }
1095                 }
1096
1097                 //calculate rcumul
1098                 for(int a = 0; a < numComp; a++) {
1099                         float rcumul = 1.0000;
1100                         //this loop fills the cumulative maps and put 0.0000 in the score freq map to make it easier to print.
1101                         for (map<double,double>::iterator it = validScores.begin(); it != validScores.end(); it++) {
1102                                 //make rscoreFreq map and rCumul
1103                                 map<double,double>::iterator it2 = rScoreFreq[a].find(it->first);
1104                                 rCumul[a][it->first] = rcumul;
1105                                 //get percentage of random trees with that info
1106                                 if (it2 != rScoreFreq[a].end()) {  rScoreFreq[a][it->first] /= iters; rcumul-= it2->second;  }
1107                                 else { rScoreFreq[a][it->first] = 0.0000; } //no random trees with that score
1108                         }
1109                 }
1110
1111         }
1112         catch(exception& e) {
1113                 m->errorOut(e, "UnifracWeightedCommand", "calculateFreqsCums");
1114                 exit(1);
1115         }
1116 }
1117 /***********************************************************/
1118
1119
1120
1121
1122