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["weighted"].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 remainingPairs = namesOfGroupCombos.size();
 703         int startIndex = 0;
 704         for (int remainingProcessors = processors; remainingProcessors > 0; remainingProcessors--) {
 705             int numPairs = remainingPairs; //case for last processor
 706             if (remainingProcessors != 1) { numPairs = ceil(remainingPairs / remainingProcessors); }
 707             lines.push_back(linePair(startIndex, numPairs)); //startIndex, numPairs
 708             startIndex = startIndex + numPairs;
 709             remainingPairs = remainingPairs - numPairs;
 710         }
 711
 712
 713
 714         //get scores for random trees
 715         for (int j = 0; j < iters; j++) {
 716 //#if defined (__APPLE__) || (__MACH__) || (linux) || (__linux) || (__linux__) || (__unix__) || (__unix)
 717             //if(processors == 1){
 718               //  driver(thisTree,  namesOfGroupCombos, 0, namesOfGroupCombos.size(),  rScores);
 719            // }else{
 720                 createProcesses(thisTree,  namesOfGroupCombos, rScores);
 721            // }
 722 //#else
 723             //driver(thisTree, namesOfGroupCombos, 0, namesOfGroupCombos.size(), rScores);
 724 //#endif
 725
 726
 727             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; }
 728
 729         }
 730         lines.clear();
 731
 732         //find the signifigance of the score for summary file
 733         for (int f = 0; f < numComp; f++) {
 734             //sort random scores
 735             sort(rScores[f].begin(), rScores[f].end());
 736
 737             //the index of the score higher than yours is returned
 738             //so if you have 1000 random trees the index returned is 100
 739             //then there are 900 trees with a score greater then you.
 740             //giving you a signifigance of 0.900
 741             int index = findIndex(usersScores[f], f);    if (index == -1) { m->mothurOut("error in UnifracWeightedCommand"); m->mothurOutEndLine(); exit(1); } //error code
 742
 743             //the signifigance is the number of trees with the users score or higher
 744             WScoreSig.push_back((iters-index)/(float)iters);
 745         }
 746
 747         //out << "Tree# " << i << endl;
 748         calculateFreqsCumuls();
 749         printWeightedFile();
 750
 751         delete output;
 752
 753         return 0;
 754     }
 755         catch(exception& e) {
 756                 m->errorOut(e, "UnifracWeightedCommand", "runRandomCalcs");
 757                 exit(1);
 758         }
 759 }
 760 /**************************************************************************************************/
 761
 762 int UnifracWeightedCommand::createProcesses(Tree* t, vector< vector<string> > namesOfGroupCombos, vector< vector<double> >& scores) {
 763         try {
 764         int process = 1;
 765                 vector<int> processIDS;
 766                 EstOutput results;
 767
 768 #if defined (__APPLE__) || (__MACH__) || (linux) || (__linux) || (__linux__) || (__unix__) || (__unix)
 769                 //loop through and create all the processes you want
 770                 while (process != processors) {
 771                         pid_t pid = fork();
 772
 773                         if (pid > 0) {
 774                                 processIDS.push_back(pid);  //create map from line number to pid so you can append files in correct order later
 775                                 process++;
 776                         }else if (pid == 0){
 777                                 driver(t, namesOfGroupCombos, lines[process].start, lines[process].num, scores);
 778
 779                                 //pass numSeqs to parent
 780                                 ofstream out;
 781                                 string tempFile = outputDir + m->mothurGetpid(process) + ".weightedcommand.results.temp";
 782                                 m->openOutputFile(tempFile, out);
 783                                 for (int i = lines[process].start; i < (lines[process].start + lines[process].num); i++) { out << scores[i][(scores[i].size()-1)] << '\t';  } out << endl;
 784                                 out.close();
 785
 786                                 exit(0);
 787                         }else {
 788                                 m->mothurOut("[ERROR]: unable to spawn the necessary processes."); m->mothurOutEndLine();
 789                                 for (int i = 0; i < processIDS.size(); i++) { kill (processIDS[i], SIGINT); }
 790                                 exit(0);
 791                         }
 792                 }
 793
 794                 driver(t, namesOfGroupCombos, lines[0].start, lines[0].num, scores);
 795
 796                 //force parent to wait until all the processes are done
 797                 for (int i=0;i<(processors-1);i++) {
 798                         int temp = processIDS[i];
 799                         wait(&temp);
 800                 }
 801
 802                 //get data created by processes
 803                 for (int i=0;i<(processors-1);i++) {
 804
 805                         ifstream in;
 806                         string s = outputDir + toString(processIDS[i]) + ".weightedcommand.results.temp";
 807                         m->openInputFile(s, in);
 808
 809                         double tempScore;
 810                         for (int j = lines[(i+1)].start; j < (lines[(i+1)].start + lines[(i+1)].num); j++) { in >> tempScore; scores[j].push_back(tempScore); }
 811                         in.close();
 812                         m->mothurRemove(s);
 813                 }
 814 #else
 815         //fill in functions
 816         vector<weightedRandomData*> pDataArray;
 817                 DWORD   dwThreadIdArray[processors-1];
 818                 HANDLE  hThreadArray[processors-1];
 819         vector<CountTable*> cts;
 820         vector<Tree*> trees;
 821
 822                 //Create processor worker threads.
 823                 for( int i=1; i<processors; i++ ){
 824             CountTable* copyCount = new CountTable();
 825             copyCount->copy(ct);
 826             Tree* copyTree = new Tree(copyCount);
 827             copyTree->getCopy(t);
 828
 829             cts.push_back(copyCount);
 830             trees.push_back(copyTree);
 831
 832             vector< vector<double> > copyScores = rScores;
 833
 834             weightedRandomData* tempweighted = new weightedRandomData(m, lines[i].start, lines[i].num, namesOfGroupCombos, copyTree, copyCount, includeRoot, copyScores);
 835                         pDataArray.push_back(tempweighted);
 836                         processIDS.push_back(i);
 837
 838                         hThreadArray[i-1] = CreateThread(NULL, 0, MyWeightedRandomThreadFunction, pDataArray[i-1], 0, &dwThreadIdArray[i-1]);
 839                 }
 840
 841                 driver(t, namesOfGroupCombos, lines[0].start, lines[0].num, scores);
 842
 843                 //Wait until all threads have terminated.
 844                 WaitForMultipleObjects(processors-1, hThreadArray, TRUE, INFINITE);
 845
 846                 //Close all thread handles and free memory allocations.
 847                 for(int i=0; i < pDataArray.size(); i++){
 848             for (int j = pDataArray[i]->start; j < (pDataArray[i]->start+pDataArray[i]->num); j++) {
 849                 scores[j].push_back(pDataArray[i]->scores[j][pDataArray[i]->scores[j].size()-1]);
 850             }
 851                         delete cts[i];
 852             delete trees[i];
 853                         CloseHandle(hThreadArray[i]);
 854                         delete pDataArray[i];
 855                 }
 856
 857
 858 #endif
 859
 860         return 0;
 861         }
 862         catch(exception& e) {
 863                 m->errorOut(e, "UnifracWeightedCommand", "createProcesses");
 864                 exit(1);
 865         }
 866 }
 867
 868 /**************************************************************************************************/
 869 int UnifracWeightedCommand::driver(Tree* t, vector< vector<string> > namesOfGroupCombos, int start, int num, vector< vector<double> >& scores) {
 870  try {
 871                 Tree* randT = new Tree(ct);
 872
 873         Weighted weighted(includeRoot);
 874
 875                 for (int h = start; h < (start+num); h++) {
 876
 877                         if (m->control_pressed) { return 0; }
 878
 879                         //initialize weighted score
 880                         string groupA = namesOfGroupCombos[h][0];
 881                         string groupB = namesOfGroupCombos[h][1];
 882
 883                         //copy T[i]'s info.
 884                         randT->getCopy(t);
 885
 886                         //create a random tree with same topology as T[i], but different labels
 887                         randT->assembleRandomUnifracTree(groupA, groupB);
 888
 889                         if (m->control_pressed) { delete randT;  return 0;  }
 890
 891                         //get wscore of random tree
 892                         EstOutput randomData = weighted.getValues(randT, groupA, groupB);
 893
 894                         if (m->control_pressed) { delete randT;  return 0;  }
 895
 896                         //save scores
 897                         scores[h].push_back(randomData[0]);
 898                 }
 899
 900                 delete randT;
 901
 902                 return 0;
 903
 904         }
 905         catch(exception& e) {
 906                 m->errorOut(e, "UnifracWeightedCommand", "driver");
 907                 exit(1);
 908         }
 909 }
 910 /***********************************************************/
 911 void UnifracWeightedCommand::printWeightedFile() {
 912         try {
 913                 vector<double> data;
 914                 vector<string> tags;
 915                 tags.push_back("Score"); tags.push_back("RandFreq"); tags.push_back("RandCumul");
 916
 917                 for(int a = 0; a < numComp; a++) {
 918                         output->initFile(groupComb[a], tags);
 919                         //print each line
 920                         for (map<double,double>::iterator it = validScores.begin(); it != validScores.end(); it++) {
 921                                 data.push_back(it->first);  data.push_back(rScoreFreq[a][it->first]); data.push_back(rCumul[a][it->first]);
 922                                 output->output(data);
 923                                 data.clear();
 924                         }
 925                         output->resetFile();
 926                 }
 927         }
 928         catch(exception& e) {
 929                 m->errorOut(e, "UnifracWeightedCommand", "printWeightedFile");
 930                 exit(1);
 931         }
 932 }
 933
 934
 935 /***********************************************************/
 936 void UnifracWeightedCommand::printWSummaryFile() {
 937         try {
 938                 //column headers
 939                 outSum << "Tree#" << '\t' << "Groups" << '\t' << "WScore" << '\t';
 940                 m->mothurOut("Tree#\tGroups\tWScore\t");
 941                 if (random) { outSum << "WSig"; m->mothurOut("WSig"); }
 942                 outSum << endl; m->mothurOutEndLine();
 943
 944                 //format output
 945                 outSum.setf(ios::fixed, ios::floatfield); outSum.setf(ios::showpoint);
 946
 947                 //print each line
 948                 int count = 0;
 949                 for (int i = 0; i < T.size(); i++) {
 950                         for (int j = 0; j < numComp; j++) {
 951                                 if (random) {
 952                                         if (WScoreSig[count] > (1/(float)iters)) {
 953                                                 outSum << setprecision(6) << i+1 << '\t' << groupComb[j] << '\t' << utreeScores[count] << '\t' << setprecision(itersString.length()) << WScoreSig[count] << endl;
 954                                                 cout << setprecision(6) << i+1 << '\t' << groupComb[j] << '\t' << utreeScores[count] << '\t' << setprecision(itersString.length()) << WScoreSig[count] << endl;
 955                                                 m->mothurOutJustToLog(toString(i+1) +"\t" + groupComb[j] +"\t" + toString(utreeScores[count]) +"\t" +  toString(WScoreSig[count]) + "\n");
 956                                         }else{
 957                                                 outSum << setprecision(6) << i+1 << '\t' << groupComb[j] << '\t' << utreeScores[count] << '\t' << setprecision(itersString.length()) << "<" << (1/float(iters)) << endl;
 958                                                 cout << setprecision(6) << i+1 << '\t' << groupComb[j] << '\t' << utreeScores[count] << '\t' << setprecision(itersString.length()) << "<" << (1/float(iters)) << endl;
 959                                                 m->mothurOutJustToLog(toString(i+1) +"\t" + groupComb[j] +"\t" + toString(utreeScores[count]) +"\t<" +  toString((1/float(iters))) + "\n");
 960                                         }
 961                                 }else{
 962                                         outSum << setprecision(6) << i+1 << '\t' << groupComb[j] << '\t' << utreeScores[count] << endl;
 963                                         cout << setprecision(6) << i+1 << '\t' << groupComb[j] << '\t' << utreeScores[count]  << endl;
 964                                         m->mothurOutJustToLog(toString(i+1) +"\t" + groupComb[j] +"\t" + toString(utreeScores[count]) +"\n");
 965                                 }
 966                                 count++;
 967                         }
 968                 }
 969                 outSum.close();
 970         }
 971         catch(exception& e) {
 972                 m->errorOut(e, "UnifracWeightedCommand", "printWSummaryFile");
 973                 exit(1);
 974         }
 975 }
 976 /***********************************************************/
 977 void UnifracWeightedCommand::createPhylipFile() {
 978         try {
 979                 int count = 0;
 980                 //for each tree
 981                 for (int i = 0; i < T.size(); i++) {
 982
 983             string phylipFileName;
 984                         map<string, string> variables;
 985             variables["[filename]"] = outputDir + m->getSimpleName(treefile);
 986             variables["[tag]"] = toString(i+1);
 987             if ((outputForm == "lt") || (outputForm == "square")) {
 988                 variables["[tag2]"] = "weighted.phylip";
 989                 phylipFileName = getOutputFileName("phylip",variables);
 990                 outputNames.push_back(phylipFileName); outputTypes["phylip"].push_back(phylipFileName);
 991             }else { //column
 992                 variables["[tag2]"] = "weighted.column";
 993                 phylipFileName = getOutputFileName("column",variables);
 994                 outputNames.push_back(phylipFileName); outputTypes["column"].push_back(phylipFileName);
 995             }
 996
 997
 998                         ofstream out;
 999                         m->openOutputFile(phylipFileName, out);
1000
1001                         if ((outputForm == "lt") || (outputForm == "square")) {
1002                                 //output numSeqs
1003                                 out << m->getNumGroups() << endl;
1004                         }
1005
1006                         //make matrix with scores in it
1007                         vector< vector<float> > dists;  dists.resize(m->getNumGroups());
1008                         for (int i = 0; i < m->getNumGroups(); i++) {
1009                                 dists[i].resize(m->getNumGroups(), 0.0);
1010                         }
1011
1012                         //flip it so you can print it
1013                         for (int r=0; r<m->getNumGroups(); r++) {
1014                                 for (int l = 0; l < r; l++) {
1015                                         dists[r][l] = utreeScores[count];
1016                                         dists[l][r] = utreeScores[count];
1017                                         count++;
1018                                 }
1019                         }
1020
1021                         //output to file
1022                         for (int r=0; r<m->getNumGroups(); r++) {
1023                                 //output name
1024                                 string name = (m->getGroups())[r];
1025                                 if (name.length() < 10) { //pad with spaces to make compatible
1026                                         while (name.length() < 10) {  name += " ";  }
1027                                 }
1028
1029                                 if (outputForm == "lt") {
1030                                         out << name << '\t';
1031
1032                                         //output distances
1033                                         for (int l = 0; l < r; l++) {   out  << dists[r][l] << '\t';  }
1034                                         out << endl;
1035                                 }else if (outputForm == "square") {
1036                                         out << name << '\t';
1037
1038                                         //output distances
1039                                         for (int l = 0; l < m->getNumGroups(); l++) {   out  << dists[r][l] << '\t';  }
1040                                         out << endl;
1041                                 }else{
1042                                         //output distances
1043                                         for (int l = 0; l < r; l++) {
1044                                                 string otherName = (m->getGroups())[l];
1045                                                 if (otherName.length() < 10) { //pad with spaces to make compatible
1046                                                         while (otherName.length() < 10) {  otherName += " ";  }
1047                                                 }
1048
1049                                                 out  << name << '\t' << otherName << dists[r][l] << endl;
1050                                         }
1051                                 }
1052                         }
1053                         out.close();
1054                 }
1055         }
1056         catch(exception& e) {
1057                 m->errorOut(e, "UnifracWeightedCommand", "createPhylipFile");
1058                 exit(1);
1059         }
1060 }
1061 /***********************************************************/
1062 int UnifracWeightedCommand::findIndex(float score, int index) {
1063         try{
1064                 for (int i = 0; i < rScores[index].size(); i++) {
1065                         if (rScores[index][i] >= score) {       return i;       }
1066                 }
1067                 return rScores[index].size();
1068         }
1069         catch(exception& e) {
1070                 m->errorOut(e, "UnifracWeightedCommand", "findIndex");
1071                 exit(1);
1072         }
1073 }
1074
1075 /***********************************************************/
1076
1077 void UnifracWeightedCommand::calculateFreqsCumuls() {
1078         try {
1079                 //clear out old tree values
1080                 rScoreFreq.clear();
1081                 rScoreFreq.resize(numComp);
1082                 rCumul.clear();
1083                 rCumul.resize(numComp);
1084                 validScores.clear();
1085
1086                 //calculate frequency
1087                 for (int f = 0; f < numComp; f++) {
1088                         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...
1089                                 validScores[rScores[f][i]] = rScores[f][i];
1090                                 map<double,double>::iterator it = rScoreFreq[f].find(rScores[f][i]);
1091                                 if (it != rScoreFreq[f].end()) {
1092                                         rScoreFreq[f][rScores[f][i]]++;
1093                                 }else{
1094                                         rScoreFreq[f][rScores[f][i]] = 1;
1095                                 }
1096                         }
1097                 }
1098
1099                 //calculate rcumul
1100                 for(int a = 0; a < numComp; a++) {
1101                         float rcumul = 1.0000;
1102                         //this loop fills the cumulative maps and put 0.0000 in the score freq map to make it easier to print.
1103                         for (map<double,double>::iterator it = validScores.begin(); it != validScores.end(); it++) {
1104                                 //make rscoreFreq map and rCumul
1105                                 map<double,double>::iterator it2 = rScoreFreq[a].find(it->first);
1106                                 rCumul[a][it->first] = rcumul;
1107                                 //get percentage of random trees with that info
1108                                 if (it2 != rScoreFreq[a].end()) {  rScoreFreq[a][it->first] /= iters; rcumul-= it2->second;  }
1109                                 else { rScoreFreq[a][it->first] = 0.0000; } //no random trees with that score
1110                         }
1111                 }
1112
1113         }
1114         catch(exception& e) {
1115                 m->errorOut(e, "UnifracWeightedCommand", "calculateFreqsCums");
1116                 exit(1);
1117         }
1118 }
1119 /***********************************************************/
1120
1121
1122
1123
1124