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